# Méthode généralisée de partage des poids et intégration de données pour l'amélioration de la précision des estimateurs de trafic postal en France

Estelle Medous

► To cite this version:

Estelle Medous. Méthode généralisée de partage des poids et intégration de données pour l'amélioration de la précision des estimateurs de trafic postal en France. Méthodologie [stat.ME]. Université Toulouse Capitole, 2023. Français. NNT: . tel-04417405

HAL Id: tel-04417405

https://theses.hal.science/tel-04417405

Submitted on 25 Jan 2024

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

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

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

Université
de Toulouse

THÈSE

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

En vue de l'obtention du

DOCTORAT DE L'UNIVERSITÉ DE TOULOUSE

Délivré par : l'Université Toulouse Capitole (UT Capitole)

Présentée et soutenue le 14/12/2023 par :
Estelle Medous

Méthode généralisée de partage des poids et intégration de données pour l'amélioration de la précision des estimateurs de trafic postal en France

JURY

<table><tr><td>PATRICE BERTAIL</td><td>PR, Univ. Paris Nanterre</td><td>Rapporteur</td></tr><tr><td>YVES TILLÉ</td><td>PR, Univ. de Neuchâtel</td><td>Rapporteur</td></tr><tr><td>PIERRE-ANDRE CORNILLON</td><td>MCF, Univ. Rennes 2</td><td>Examinateur</td></tr><tr><td>ERIC GAUTIER</td><td>PR, Univ. Toulouse Capitole</td><td>Président du jury</td></tr><tr><td>ANNE RUIZ-GAZEN</td><td>PR, Univ. Toulouse Capitole</td><td>Directrice</td></tr><tr><td>CAMELIA GOGA</td><td>PR, Univ. de Franche-Comté</td><td>Co-Directrice</td></tr><tr><td>JEAN-FRANÇOIS BEAUMONT</td><td>Chef Recherche, Statistique Canada</td><td>Invité</td></tr><tr><td>ALAIN DESSERTAINE</td><td>Expert Stratégie DATA, La Poste</td><td>Invité</td></tr><tr><td>PAULINE PUECH</td><td>Resp. Statistique, La Poste</td><td>Invitée</td></tr></table>

École doctorale et spécialité :

MITT : Domaine Mathématiques : Mathématiques appliquées

Unité de Recherche :

Toulouse School of Economics (TSE-R)

Directrices de Thèse :

Anne RUIZ-GAZEN et Camelia GOGA

Rapporteurs :

Patrice BERTAIL et Yves TILLÉ
Remerciements

Je tiens en premier lieu à remercier mes directrices de thèse, Anne Ruiz-Gazen et Camelia Goga, pour m'avoir accompagnée tout au long de ces quatre ans et m'avoir transmis une partie de leur connaissance en théorie des sondages ainsi qu'en rédaction, bien que je sois encore loin de leur niveau dans ces deux domaines. Je voudrais aussi remercier mes encadrants à la Poste, Alain Dessertaine et Pauline Puech, qui m'ont formée et encouragée tout au long de mon stage de Master puis pendant la thèse. Grace à eux j'ai pu découvrir le monde des sondages en entreprise et mettre en pratique la théorie.

Je souhaite aussi remercier mes rapporteurs de thèse, Partice Bertail et Yves Tillé, pour avoir pris le temps de lire mon manuscrit et pour les précieux retours qu'ils m'ont apportés. Je tiens aussi à remercier Pierre-Andre Cornillon et Eric Gautier pour avoir accepté de faire partie de mon jury de thèse.

Je remercie aussi Jean-François Beaumont pour m'avoir fait découvrir Ottawa et pour ses remarques et encouragements.

Je remercie aussi La Poste, qui a financé cette thèse et le stage de Master la précédant. Je remercie particulièrement le service Statistique de La Poste pour m'avoir accueillie et encadrée dans le monde entreprenarial malgré les difficultés rencontrées, et les membres de l'UMR TSE-R pour leur accueil et leur échanges.

Je voudrais remercier ma famille pour leur soutien. Je remercie particulièrement mon frère Charles et un ami de longue date, Guillaume, qui ont fait les mêmes choix d'orientation que moi, pour leurs discussions édifiantes. Je remercie aussi Pepper pour avoir daigné m'ac-compagner un peu plus longtemps et je remercie Amalric pour avoir été là dans les moments les plus difficiles.
<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Pepper.</figcaption></figure>
# Table des matières

<table><tr><td>Introduction</td><td>9</td></tr><tr><td>Acronymes</td><td>13</td></tr><tr><td>1 Rappels et contexte</td><td>15</td></tr><tr><td>1.1 Théorie des sondages</td><td>15</td></tr><tr><td>1.1.1 Principes</td><td>15</td></tr><tr><td>1.1.2 Échantillonnage</td><td>17</td></tr><tr><td>1.1.3 Estimation</td><td>24</td></tr><tr><td>1.1.4 Sondage indirect</td><td>33</td></tr><tr><td>1.2 La Poste</td><td>37</td></tr><tr><td>1.2.1 Plans de sondage SYCI 1 et 2</td><td>38</td></tr><tr><td>1.2.2 Données auxiliaires</td><td>44</td></tr><tr><td>1.3 Contributions</td><td>50</td></tr><tr><td>1.3.1 Sondage indirect et MGPP</td><td>51</td></tr><tr><td>1.3.2 Intégration statistique de données</td><td>62</td></tr><tr><td>Sondage indirect et MGPP</td><td>69</td></tr><tr><td>2 Indirect sampling</td><td>69</td></tr><tr><td>2.1 Introduction</td><td>70</td></tr><tr><td>2.2 Indirect sampling</td><td>72</td></tr><tr><td>2.2.1 GWSM</td><td>72</td></tr><tr><td>2.2.2 MtO links and optimal weight links</td><td>75</td></tr><tr><td>2.2.3 Comparison of direct and optimal indirect MtO sampling designs</td><td>78</td></tr><tr><td>2.3 Double indirect sampling</td><td>79</td></tr><tr><td>2.3.1 Double GWSM</td><td>79</td></tr><tr><td>2.3.2 MtO links</td><td>81</td></tr></table>
<table><tbody><tr><td>2.3.3</td><td>MtO-MtO links and double standardization</td><td>81</td></tr><tr><td>2.4</td><td>Simulation study</td><td>84</td></tr><tr><td>2.4.1</td><td>Population and link setups</td><td>84</td></tr><tr><td>2.4.2</td><td>Sampling designs and GWSM estimators</td><td>86</td></tr><tr><td>2.4.3</td><td>Measures of efficiency and results</td><td>88</td></tr><tr><td>2.5</td><td>Application to the French post data</td><td>90</td></tr><tr><td>2.6</td><td>Conclusion</td><td>92</td></tr><tr><td>2.7</td><td>Appendix</td><td>93</td></tr><tr><td><b>3</b></td><td><b>Optimal Weights for double GWSM</b></td><td><b>99</b></td></tr><tr><td>3.1</td><td>Introduction</td><td>100</td></tr><tr><td>3.2</td><td>Double Generalized Weight Share Method</td><td>102</td></tr><tr><td>3.2.1</td><td>Indirect sampling and simple GWSM estimator</td><td>102</td></tr><tr><td>3.2.2</td><td>Double GWSM estimator</td><td>104</td></tr><tr><td>3.3</td><td>Optimality of double GWSM</td><td>106</td></tr><tr><td>3.4</td><td>Alternative to double GWSM optimal link weights</td><td>113</td></tr><tr><td>3.4.1</td><td>Alternative link weights using auxiliary information</td><td>113</td></tr><tr><td>3.4.2</td><td>How to choose the best auxiliary variable</td><td>114</td></tr><tr><td>3.5</td><td>Simulations</td><td>116</td></tr><tr><td>3.5.1</td><td>Populations</td><td>116</td></tr><tr><td>3.5.2</td><td>Parameters and results</td><td>117</td></tr><tr><td>3.5.3</td><td>Perspectives for La Poste application</td><td>120</td></tr><tr><td>3.6</td><td>Conclusion</td><td>121</td></tr><tr><td>3.7</td><td>Appendix</td><td>121</td></tr></tbody></table>

Intégration statistique de données 135

4 QR prediction for SDI 135

<table><tbody><tr><td>4.1</td><td>Introduction</td><td>136</td></tr><tr><td>4.2</td><td>Study variable observed in both samples</td><td>137</td></tr><tr><td>4.3</td><td>study variable unobserved in NP sample</td><td>140</td></tr><tr><td>4.4</td><td>Asymptotic properties and variance estimation of QR predictors</td><td>144</td></tr><tr><td>4.4.1</td><td>Bias properties of T̂<sup>(Qπ)</sup></td><td>145</td></tr><tr><td>4.4.2</td><td>Asymptotic variance and variance estimation of T̂<sup>(Qπ)</sup></td><td>145</td></tr><tr><td>4.5</td><td>Simulations</td><td>147</td></tr><tr><td>4.5.1</td><td>Populations and setups</td><td>147</td></tr></tbody></table>
<table><tr><td>4.5.2</td><td>Results</td><td>149</td></tr><tr><td>4.6</td><td>Application to La Poste data</td><td>155</td></tr><tr><td>4.6.1</td><td>Data presentation</td><td>155</td></tr><tr><td>4.6.2</td><td>Results</td><td>156</td></tr><tr><td>4.7</td><td>Conclusion</td><td>158</td></tr><tr><td>4.8</td><td>Supplementary material</td><td>161</td></tr><tr><td>4.8.1</td><td>Assumptions</td><td>161</td></tr><tr><td>4.8.2</td><td>Development of the results presented in section 4.4.2</td><td>161</td></tr><tr><td colspan="2">Conclusions et perspectives</td><td>165</td></tr><tr><td colspan="2">Bibliographie</td><td>167</td></tr></table>
# Introduction

La Poste possède un observatoire du trafic postal français qui lui permet, entre autres, de s'adapter à la réalité du terrain et de développer de nouveaux produits. Cet observatoire utilise une étude par sondage nommée **SYCI 2** afin de recueillir des données comme le volume et la nature du trafic postal, les délais d'acheminement, ou encore les flux géographiques et économiques.

Actuellement, La Poste cherche à réduire la taille des échantillons prélevés pour **SYCI 2**, ce qui compromet la précision des estimateurs du trafic total. Mon objectif de thèse consiste à comprendre et améliorer la méthodologie de **SYCI 2** pour maintenir la précision des estimateurs de La Poste malgré la réduction des tailles d'échantillons. Mes travaux de recherche s'organisent en deux axes.

Le premier axe consiste à améliorer la méthodologie actuelle, basée sur l'utilisation d'une **Méthode Généralisée de Partage des Poids**. Cette méthode, introduite par Deville and Lavallée (2006) puis développée dans Lavallée (2007), permet, sous certaines contraintes, d'obtenir un estimateur du total, dit optimal, qui est le plus précis possible en terme de variance, quelle que soit la variable d'intérêt. Bien que de nombreux travaux existent sur ce sujet, tels que Deville and Maumy-Bertrand (2006), De Vitiis et al. (2014), Kiesl (2016) et Haziza and Beaumont (2017), la question de la précision de la **MGPP** n'est abordée que dans Deville and Lavallée (2006) et Lavallée (2007), où les auteurs montrent l'existence d'estimateurs **MGPP** optimaux dans quelques cas particuliers. Mes recherches visent à identifier des situations pour lesquelles un estimateur **MGPP** optimal existe et à proposer des alternatives applicables par La Poste en utilisant les données auxiliaires récoltées dans le cadre de l'étude **SYCI 2**.

Le deuxième axe porte sur l'utilisation des bases de données massives dont dispose La Poste. Certaines informations sur les objets postaux (lettres, colis ...) triés par machines de tri sont automatiquement récupérées. La Poste souhaite développer une nouvelle méthodologie pour **SYCI 2** qui prenne en compte ces données par **Intégration Statistique de Données (ISD)**. L' **ISD** est un sujet récent mais bien étudié dans la littérature, comme en témoignent des études telles que celles menées par Kim and Tam (2021), Kim et al. (2021) ainsi que les revues
exhaustives des différentes méthodes d' ISD réalisées par Beaumont (2020), Rao (2021), Kim (2022), Wu (2022) et Yang and Kim (2020). L' ISD consiste à utiliser, d'une part, une base de données non probabiliste et d'autre part un échantillon probabiliste pour obtenir des estimations de paramètres sur la population d'intérêt, ici les objets postaux. Cependant, les méthodes étudiées dans la littérature requièrent que les valeurs des variables d'intérêt soient connues pour tout individu contenu dans la base non probabiliste. Les données collectées par machine de tri ne remplissent pas cette condition, car certaines variables d'intérêt, comme le poids et les dimensions d'un objet, ne sont pas récoltées par les machines. La Poste ne peut donc pas utiliser ces méthodes. Le but de mes recherches, dans ce second axe, est de proposer de nouvelles méthodes d' ISD adaptées aux données de La Poste.

Le chapitre 1 présente en détail mon sujet de thèse et la problématique de La Poste. Je commence par rappeler les bases de la théorie des sondages, nécessaires à la compréhension de mes travaux. Une présentation détaillée de l'étude SYCI 2 et des données de La Poste permet de mieux saisir les enjeux pour La Poste. Finalement, je résume les principales contributions apportées dans ma thèse.

Les chapitres 2 et 3 détaillent mes travaux sur le premier axe de recherche. Dans le chapitre 2, on trouve l'article Medous et al. (2023a), publié dans *The Annals of Applied Statistics*. On définit des conditions d'existence pour l'estimateur MGPP optimal puis on introduit le concept de MGPP double, qui est le cas particulier de MGPP (aussi appelée MGPP simple) utilisé par La Poste. Ce cas particulier offre des avantages en termes de faisabilité qui sont étudiés dans le chapitre 2.

La MGPP double mène cependant à une perte de précision comparée à la MGPP simple, comme l'a expérimenté La Poste. L'existence d'estimateurs optimaux pour la MGPP double est discutée dans le chapitre 3, qui fera l'objet d'un article que je compte soumettre prochainement en tant que seule autrice. Je propose aussi dans ce travail des estimateurs alternatifs pour les cas où l'estimateur optimal ne peut être utilisé.

Le chapitre 4 résume mes travaux sur le deuxième axe de recherche. Il contient l'article Medous et al. (2023b), accepté dans *Survey Methodology/Techniques d'Enquête* et qui porte sur l'intégration statistique de données. On commence par étudier certaines méthodes de la littérature pour ensuite proposer des méthodes adaptées aux données observées à La Poste.

Les différents travaux présentés dans ce mémoire ont donné lieu à des publications :

- E. Medous, C. Goga, A. Ruiz-Gazen, J.F. Beaumont, A. Dessertaine and P.Puech (2023), Many-to-One indirect sampling with application to the French postal traffic estimation,

*The Annals of Applied Statistics*, Institute of Mathematical Statistics, volume 17(1) pages 838–859,
- E. Medous, C. Goga, A. Ruiz-Gazen, J.F. Beaumont, A. Dessertaine and P.Puech (2023), QR Prediction for Statistical Data Integration,

*Survey Methodology/Techniques d'Enquête*, Statistique Canada, Numéro consacré au 11e Colloque International Francophone sur les Sondages, à paraître,

et je les ai présentés lors des conférences suivantes :

- Comparaison des sondages indirects simple et double. Application à l'estimation du trafic postal en France.

- Computational and Methodological Statistics (CMStat), en ligne, 19-21 Decembre 2020,

- *Journées de Statistique, Nice (en ligne)*, 7-11 Juin 2021.

- Pros and Cons of the Double Indirect Sampling for “Many to One” links.

*Conference in Honour of Fred Smith & Chris Skinner, Southampton, Angleterre (en ligne)*, 7-9 Juillet 2021.

- Estimation du trafic postal en France : perte de précision dûe à l'utilisation d'un sondage indirect double.

*Colloque francophone sur les sondages, Bruxelles, Belgique, 6-8 Octobre 2021.*

- Une approche par prédiction pour l'intégration de données.

*Forum des Jeunes Mathématicien.nes, Besançon, 8-10 Décembre 2021.*

- Introduction à l'intégration de données en sondages.

*Rencontre des Jeunes Statisticiens, Porquerolles, 3-7 Avril 2022.*

- Improving finite population inference by data integration & Statistical data integration using a prediction approach.

*5th International Conference on Econometrics and Statistics, Kyoto, Japon (en ligne)*, 4-6 Juin 2022 (avec A. Ruiz-Gazen).

- Optimality of the double Generalized Weight Share Method and alternatives.

- *Colloque francophone sur les sondages, Aubervilliers, 21-24 mars 2023,*

- *Journées de Statistiques, Bruxelles, Belgique, 3-7 Juillet 2023.*

ainsi qu'aux séminaires suivants :

- Comparaison des sondages indirects simple et double. Application à l'estimation du trafic postal en France.

*Doctoral Workshop on Decision Mathematics and Statistics, Toulouse, 6 Juillet 2021.*

- Pros and Cons of the Double Indirect Sampling for “Many to One” links.

*Seminar of the University of Jyväskylä, Jyväskylä, Finlande (en ligne)*, 5 Novembre 2021.

- Many-to-One indirect sampling with application to the French postal traffic estimation.
Séminaire MAD-Stat, Toulouse, 10 Mars 2022, (avec A. Ruiz-Gazen).

- A prediction approach for statistical data integration in survey sampling.

Doctoral Workshop on Decision Mathematics and Statistics, Toulouse, 16 Juin 2022.

- Optimality of the double Generalized Weight Share Method and alternatives.

Séminaire du Laboratoire de mathématiques de Besançon, 31 Janvier 2023.

Mes travaux portant sur le sondage indirect ont été récompensés par un prix lors de la
*Student competition in association with the Conference in Honour of Fred Smith & Chris*
*Skinner*, du 7 au 9 Juillet 2021.
# Acronymes

**DGWSM** Double Generalized Weight Share Method.

**EQM** Erreur Quadratique Moyenne.

**GWSM** Generalized Weight Share Method.

**HT** Horvitz-Thompson.

**ISD** Intégration Statistique de Données.

**MGPP** Méthode Généralisée de Partage des Poids.

**MtO** Many-to-One.

**MtO-MtO** Many-to-One-Many-to-One.

**RAO** Référentiel des Adresses Organisées.

**RRMSE** Relative Root Mean Square Error.

**SASSR** Sondage Aléatoire Simple Sans Remise.

**SDI** Statistical Data Integration.

**SOJ** SOrties-Jours.

**SRSWOR** Simple Random Sampling Without Replacement.

**SSRSWOR** Stratified Simple Random Sampling Without Replacement.

**STASSR** sondage STratifié Aléatoire Simple Sans Remise.

**SYCI** SYstème de Collecte de l'Information.

**TAE** Traitement Automatisé de l'Enveloppe.

**TpU** Tous-pour-Un.

**TpU-TpU** Tous-pour-Un-Tous-pour-Un.

**UP** Unités Primaires.

**US** Unités Secondaires.
# Chapitre 1

# Rappels de théorie des sondages,
contexte et contributions

## 1.1 Théorie des sondages

Dans cette section, on introduit les notions de théorie des sondages nécessaires à la compréhension de mes travaux de thèse. On rappelle les méthodes classiques d’échantillonnage et d’estimation, ainsi que le concept d’asymptotique en théorie des sondages. On présente ensuite le sondage indirect et la méthode généralisée de partage des poids (MGPP) sur lesquels porte une grande partie de mes travaux. Pour plus de détails sur cette section, on peut se référer aux ouvrages suivants : Hájek (1981), Särndal et al. (1989), Ardilly (2006), Lavallée (2007) et Tillé (2019).

### 1.1.1 Principes

On appelle population finie de taille *N*, notée *U*, une collection de *N* individus possédant un ensemble de caractéristiques appelées variables. On peut considérer par exemple la population des chats français comme dans la figure 1.1 ou, plus communément, la population des personnes habitant en France. Ces personnes forment une population finie de taille 68 042 591 (début 2023) et sont caractérisés par leur sexe, taille, âge, adresse...

On dit qu'une population *U* est cible si on souhaite connaître certains paramètres pour cette population, comme le total ou la moyenne d'une certaine variable. On peut aussi s'intéresser à des paramètres plus complexes, comme la médiane d'une certaine variable ou les covariances entre différentes variables. Le paramètre étudié et les variables qui permettent de le calculer sont dits d'intérêt. Par exemple, on peut s'intéresser au pourcentage de femmes vivant en France, à leur revenu médian ou au nombre total de chats de compagnie en France.
<figure><img src="image_4.png" /><figcaption>FIGURE 1.1 – Population des (chats) français.</figcaption></figure>

La variable d'intérêt est respectivement le sexe d'un individu, son revenu et le nombre de chats résidant chez lui. Les variables qui ne sont pas d'intérêt, comme le niveau d'étude ou le type de logement, sont dites auxiliaires.

Soit $\theta_y$ un paramètre d'intérêt qui est une fonction d'une variable d'intérêt $y$ sur $U$. On peut calculer la valeur de ce paramètre si on observe tous les individus de la population. Cette méthode, appelée recensement, est coûteuse et complexe à mettre en place. Par exemple, pour recenser les personnes habitant en France, il faut observer chacune des 68 042 591 personnes composant la population et leur poser des questions, soit par téléphone, par mail, en face à face... Cela nécessite donc une infrastructure conséquente.

Une autre méthode, appelée sondage, consiste à sélectionner, le plus souvent de manière aléatoire, un sous-groupe d'individus $s$, appelé échantillon, comme illustré dans la figure 1.2. Les paramètres d'intérêt sont ensuite estimés à partir des valeurs des variables d'intérêt observées pour les individus échantillonnés. On appelle méthode d'échantillonnage la façon dont est tiré l'échantillon. Différentes méthodes d'échantillonnage seront présentées en section 1.1.2.

Mettre en place un sondage comporte de nombreuses étapes. Il faut d'abord identifier les objectifs de l'étude ainsi que la population et les paramètres d'intérêt. Ensuite, il faut identifier les informations auxiliaires disponibles au niveau de la population d'intérêt. Ces données peuvent servir à décider la façon dont les individus seront sélectionnés dans l'échantillon.

Les données récoltées au niveau de l'échantillon, parfois couplées avec les informations auxiliaires disponibles au niveau de la population, vont permettre d'estimer les paramètres d'intérêt. L'étape qui consiste à calculer ces estimateurs est appelée estimation et sera détaillée en section 1.1.3.

Il existe d'autres étapes dans la mise en pratique d'un sondage, comme la rédaction du questionnaire ou l'encodage des résultats. Cependant, mes travaux portant uniquement sur
<figure><img src="image_5.png" /><figcaption>FIGURE 1.2 – Échantillon de (chats) français.</figcaption></figure>

les étapes d'échantillonnage et d'estimation, je ne m'attarderai pas sur les autres étapes dans
ce chapitre.

1.1.2 Échantillonnage

Dans cette section, ainsi que dans le reste du chapitre, on considère les variables d’intérêt comme non aléatoires.

Tirer un échantillon de manière aléatoire nécessite d'identifier de manière unique tous les individus de la population. On appelle base de sondage la liste contenant tous les individus de U identifiés de manière unique, comme illustré dans la figure 1.3. La base de sondage contient généralement des informations auxiliaires sur les individus de la population, comme le lieu de naissance, le nom...

Soit S l'ensemble des échantillons possibles et $\mathcal{B}(S)$ la tribu borélienne engendrée par S. La taille d'un échantillon $s \in S$, notée $n_s$, peut être fixe ou variable. On appelle plan de sondage $p(\cdot)$ une loi de probabilités sur $(S, \mathcal{B}(S))$ telle que $p(s)$ est la probabilité d'échantillonner $s \in S$. Le plan de sondage consiste donc à attribuer une probabilité de tirage à chaque échantillon $s \in S$. Un plan de sondage est déterminé avant le tirage en utilisant les données accessibles dans la base de sondage. Dans cette section, je présente les plans de sondage classiques utilisés dans mes travaux.

Pour un individu $k \in U$, on note $I_k$ l'indicatrice d'inclusion dans un échantillon de l'individu $k \in U$, où $I_k = 1$ si k est sélectionné et 0 sinon. Les indicatrices d'inclusion sont aléatoires et un échantillon est composé de tous les individus $k \in U$ dont la valeur de $I_k$ vaut 1.

On appelle probabilité d'inclusion de premier ordre de l'individu $k \in U$, notée $\pi_k$, la probabilité que l'individu $k$ soit sélectionné dans un échantillon et probabilité d'inclusion
<figure><img src="image_6.png" /><figcaption>FIGURE 1.3 – Base de sondage des (chats) français.</figcaption></figure>

de second ordre du couple $k, k' \in U$, notée $\pi_{kk'}$, la probabilité que $k$ et $k'$ soit sélectionnés dans un échantillon. Soit $E[\cdot]$ l'espérance sous le plan de sondage $p(\cdot)$. On a $\pi_k = E[I_k]$ et $\pi_{kk'} = E[I_k I_{k'}]$, ce qui implique que $\pi_{kk} = E[I_k I_k] = E[I_k] = \pi_k$. La covariance entre $I_k$ et $I_{k'}$ est donnée par :

$$ \mathrm{Cov}(I_k, I_{k'}) = E(I_k I_{k'}) - E(I_k)E(I_{k'}) = \pi_{kk'} - \pi_k \pi_{k'}, \quad k, k' \in U. $$

Les probabilités d'inclusion de premier et second ordre permettent de construire un estimateur de $\theta_y$ et de déterminer sa précision, ce qui sera détaillé dans la section 1.1.3. Les valeurs de $\pi_k$ et $\pi_{kk'}$, $k, k' \in U$, seront données pour chaque plan présenté dans cette section.

## Sondage Aléatoire Simple Sans Remise

Le Sondage Aléatoire Simple Sans Remise (SASSR) est un plan très simple à mettre en oeuvre qui permet de tirer des échantillons sans remise d'une taille donnée $n$. Le SASSR consiste à attribuer la même probabilité d'être tiré à chaque échantillon sans remise $s$ de taille $n$ et une probabilité nulle aux autres échantillons. Comme il existe $\binom{N}{n}$ échantillons sans remise de taille $n$ dans $U$, la probabilité de tirage d'un échantillon $s$ est $p(s) = 1/\binom{N}{n}$. On peut alors calculer la probabilité que l'individu $k$ soit sélectionné dans un échantillon, qui est le ratio entre le nombre d'échantillons contenant $k$ et le nombre total d'échantillons
possibles :

$$
\pi_k = \frac{\binom{N-1}{n-1}}{\binom{N}{n}} = \frac{n}{N}, k \in U.
$$

On peut dériver de manière similaire la probabilité d'inclusion de second ordre :

$$
\pi_{kk'} = \frac{\binom{N-2}{n-2}}{\binom{N}{n}} = \frac{n(n-1)}{N(N-1)}, k \neq k' \in U.
$$

Avec un SASSR, les probabilités d'inclusion de premier ordre sont identiques pour tous les individus $k \in U$. On dit que c'est un plan de sondage à probabilités égales.

<figure><img src="image_7.png" /><figcaption>FIGURE 1.4 – Echantillons tirés via SASSR.</figcaption></figure>

Par exemple, dans la figure 1.4, on a une population composée de 7 chats français et on veut tirer un échantillon sans remise de taille 2 dans cette population. Il existe 21 échantillons de taille 2, donc la probabilité de tirer un de ces échantillons est 1/21. De là, on calcule la probabilité d'observer un chat français dans cette population, qui est de 2/7 pour tous les chats.

Plan de Poisson

Le plan de Poisson permet de tirer des échantillons sans remise tels que les indicatrices d'inclusion soient indépendantes.

On fixe la probabilité d'inclusion du premier ordre $\pi_k$ pour chaque individu $k \in U$ grâce aux informations disponibles dans la base de sondage et on sélectionne de façon indépendante et sans remise les individus avec les probabilités prédéfinies $\pi_1, \dots, \pi_N$. Si les probabilités $\pi_k, k \in U$ sont choisies fixes égales à $\pi$, le plan est dit de Bernoulli.

Pour un plan de Poisson, les indicatrices d'inclusion $I_k$ pour $k \in U$ suivent une loi de Bernoulli de paramètre $\pi_k$ et elles sont indépendantes. La taille $n_s$ de l'échantillon est
donc aléatoire. Si les probabilités $\pi_k$ sont égales à $\pi$, alors la taille de l'échantillon suit une binomiale de paramètres $N$ et $\pi$.

L'indépendance des indicatrices d'inclusion implique que pour tout $k, k' \neq k \in U$, $\pi_{kk'} = \pi_k \pi_{k'}$ et que la covariance entre deux indicatrices $I_k$ et $I_{k'}$ avec $k \neq k'$ est nulle.

<figure><img src="image_8.png" /><figcaption>FIGURE 1.5 – Échantillons tirés via plan de Poisson.</figcaption></figure>

Dans la figure 1.5, on fixe les probabilités $\pi_k$ pour chaque chat français. Trois échantillons possibles avec un plan de Poisson sont présentés en rouge. L'un est de taille 1, le deuxième de taille 2 et le troisième de taille 3.

## Plan stratifié

Si la base de sondage contient suffisamment d'information sur les individus de $U$, il est possible de partitionner $U$ dans des sous-populations disjointes $U_h$, $h = 1, \dots, H$ avec $U = \cup_{h=1}^H U_h$. Ces sous-populations s'appellent des strates. Un plan de sondage stratifié consiste à tirer de manière indépendante un échantillon $s_h$ dans chaque strate $U_h$, $h = 1, \dots, H$. L'échantillon final est donné par $s = \cup_{h=1}^H s_h$. Si le tirage au sein de chaque strate s'effectue avec un SASSR, alors le plan est appelé **sondage STatifié Aléatoire Simple Sans Remise (STASSR)**.

Ce type d'échantillonnage est utile si on souhaite par exemple effectuer des estimations pour certaines sous-populations correspondant aux strates. Par exemple, si on souhaite connaître le pourcentage de propriétaires de chat dans chaque région française. Il permet aussi d'améliorer dans certaines situations la précision des estimateurs de $\theta_y$, ce qui sera détaillé en section 1.1.3.

Soit l'ensemble de strates $U_h$, $h = 1, \dots, H$ de taille respective $N_h$ correspondant à la partition de $U$. Pour un STASSR où l'échantillon $s_h$ dans la strate $U_h$, $h = 1 \dots H$, est de
<figure><img src="image_9.png" /><figcaption>FIGURE 1.6 – Échantillon tiré via STASSR.</figcaption></figure>

taille $n_h$, la probabilité d'inclusion de premier ordre pour un individu $k \in U$ est :

$$ \pi_k = \frac{n_h}{N_h} \text{ si } k \in U_h. $$

La probabilité d'inclusion de second ordre pour deux individus $k, k' \neq k \in U$ est :

$$ \pi_{kk'} = \begin{cases} \frac{n_h(n_h-1)}{N_h(N_h-1)} & \text{si } k, k' \in U_h, \\ \frac{n_h}{N_h} \frac{n_{h'}}{N_{h'}} & \text{si } k \in U_h \text{ et } k' \in U_{h'}, h \neq h'. \end{cases} $$

Dans la figure 1.6, la population de 7 chats français est décomposée en trois strates $U_1, U_2$ et $U_3$ de taille respective 2, 3 et 2. Dans chaque strate on tire un échantillon de taille 1 par SASSR. Les probabilités d'inclusion de premier ordre sont donc 1/2 pour les individus appartenant aux strates 1 et 3, et 1/3 pour ceux appartenant à la strate 2.

## Plan à deux degrés

Le plan SASSR, de Poisson ou les plans stratifiés, nécessitent d'avoir une base de sondage des individus. Quand cette base de sondage est difficile à obtenir, on peut faire appel aux plans à plusieurs degrés. Nous allons dans la suite décrire le plan à deux degrés.

La population $U$ est découpée en $N_I$ sous groupes $U_i$, $i \in \{1, \dots, N_I\}$ de taille respective $N_i$, appelés **Unités Primaires (UP)**. Les individus $k \in U_i$, $i \in \{1, \dots, N_I\}$ sont appelés **Unités Secondaires (US)**.

Le plan à deux degrés consiste à tirer un échantillon $s_I$ d'UP avec les informations disponibles dans la base de sondage des UP. Pour chaque UP $U_i$ contenue dans $s_I$, on récolte suffisamment d'information pour former une base de sondage adéquate pour $U_i$. Cette base de sondage permet de tirer un échantillon d'US $s_i$ au sein de l'UP $U_i$. L'échantillon final $s$
est composé de tous les individus contenus dans les échantillons $s_i$, $i \in s_I$, $s = \cup_{i \in s_I} s_i$.

Pour tout $i, i' \in \{1, \dots, N_I\}$, on note $\pi_{Ii}$ la probabilité d'inclusion de l'UP $U_i$ et $\pi_{Iii'}$ la probabilité d'inclusion des UP $U_i$ et $U_{i'}$ dans un échantillon d'UP. Soient $k$ et $k'$ deux US contenues respectivement dans les UP $U_i$ et $U_{i'}$. On note $\pi_{IIk|i}$ et $\pi_{IIkk'|i,i'}$ les probabilités conditionnelles d'inclusion dans $s$ de, respectivement, l'US $k \in U_i$ et les deux US $k \in U_i$, $k' \in U_{i'}$ sachant que les UP $U_i$ et $U_{i'}$ appartiennent à l'échantillon $s_I$. Si $i = i'$, on note $\pi_{IIkk'|i,i} = \pi_{IIkk'|i}$.

On note $S_I$ l'ensemble des échantillons d'UP possibles et $p(s_I)$ la probabilité de tirer l'échantillon $s_I$. La probabilité d'inclusion dans l'échantillon final $s$ de premier ordre pour un individu $k \in U$ est égale à :

$$ \pi_k = \sum_{\substack{s_I \in S_I \\ U_i \in s_I}} p(s_I) \pi_{IIk|i} \text{ si } k \in U_i. $$

La probabilité d'inclusion dans l'échantillon final $s$ de second ordre pour deux individus $k, k' \in U$ est égale à :

$$ \pi_{kk'} = \sum_{\substack{s_I \in S_I \\ U_i, U_{i'} \in s_I}} p(s_I) \pi_{IIkk'|i,i'} \text{ si } k \in U_i, k' \in U_{i'}. $$

Pour un plan à deux degrés défini de façon générale comme ci-dessus, le calcul des probabilités d'inclusion des individus $k \in U$ peut s'avérer assez complexe. C'est pour cette raison que l'on suppose en général les conditions supplémentaires suivantes :

- le tirage des unités secondaires est effectué indépendamment au sein de chaque UP (indépendance),

- le plan de sondage utilisé pour tirer les unités secondaires ne dépend pas de l'échantillon $s_I$ tiré au premier degré (invariance).

Sous ces conditions, les probabilités d'inclusion des individus $k \in U$ dans l'échantillon final $s$ se simplifient. La probabilité d'inclusion de premier ordre pour un individu $k \in U$ devient :

$$ \pi_k = \sum_{\substack{s_I \in S_I \\ U_i \in s_I}} p(s_I) \pi_{IIk|i} = \pi_{IIk|i} \sum_{\substack{s_I \in S_I \\ U_i \in s_I}} p(s_I) = \pi_{IIk|i} \pi_{Ii} \text{ si } k \in U_i, $$

car $\pi_{IIk|i}$ ne dépend pas de $s_I$ par invariance. La probabilité d'inclusion de second ordre pour deux individus $k, k' \in U$ devient :
<figure><img src="image_10.png" /><figcaption>FIGURE 1.7 – Echantillons tirés via un plan à deux degrés avec SASSR à chaque degré.</figcaption></figure>

— Si $i = i'$, alors par invariance, $\pi_{IIkk'|i}$ ne dépend pas de $s_I$ :

$$ \pi_{kk'} = \sum_{\substack{s_I \in S_I \\ U_i \in s_I}} p(s_I) \pi_{IIkk'|i} = \pi_{IIkk'|i} \sum_{\substack{s_I \in S_I \\ U_i \in s_I}} p(s_I) = \pi_{IIkk'|i} \pi_{Ii} \text{ si } k, k' \in U_i, $$

— si $i \neq i'$, alors par indépendance $\pi_{IIkk'|i,i'} = \pi_{IIk|i} \pi_{IIk'|i'}$ et par invariance, $\pi_{IIk|i}$ et $\pi_{IIk'|i'}$ ne dépendent pas de $s_I$, alors pour tout $k \in U_i$, $k' \in U_{i'}$ :

$$ \pi_{kk'} = \sum_{\substack{s_I \in S_I \\ U_i, U_{i'} \in s_I}} p(s_I) \pi_{IIk|i} \pi_{IIk'|i'} = \pi_{IIk|i} \pi_{IIk'|i'} \sum_{\substack{s_I \in S_I \\ U_i, U_{i'} \in s_I}} p(s_I) = \pi_{IIk|i} \pi_{IIk'|i'} \pi_{Iii'}. $$

Soit un plan à deux degrés avec SASSR à chaque degré, avec $n_I$ le nombre d'UP échantillonnées au premier degré et $n_i$ le nombre d'US échantillonnées au deuxième degré au sein de l'UP $U_i$, $i \in \{1, \dots, N_I\}$. Alors les probabilités d'inclusion de premier et second ordre pour les UP sont données par :

$$ \pi_{Ii} = \frac{n_I}{N_I}, \quad \pi_{Iii'} = \frac{n_I(n_I - 1)}{N_I(N_I - 1)}, \quad i, i' \neq i \in \{1, \dots, N_I\}. $$

Les probabilités conditionnelles d'inclusion de premier et second ordre pour les US sont données par :

$$ \pi_{IIk|i} = \frac{n_i}{N_i}, \quad \pi_{IIkk'|i} = \frac{n_i(n_i - 1)}{N_i(N_i - 1)}, \quad k, k' \neq k \in U_i \text{ et } \pi_{IIkk'|i,i'} = \frac{n_i}{N_i} \frac{n_{i'}}{N_{i'}} \text{ si } k \in U_i, k' \in U_{i'}, i \neq i'. $$

Considérons maintenant l'exemple de la figure 1.7. La population est découpée en trois UP de taille 2 pour $U_1$ et $U_3$ et 3 pour $U_2$. Les chats français constituent les US. Un SASSR est utilisé à chaque degré, avec une UP échantillonnée au premier degré et une US échantillonnée
par **UP** de $s_I$. La probabilité d'observer une **UP** est donc de $1/3$ pour toutes les **UP**. Dans notre exemple, l'**UP** $U_2$ est échantillonnée. Au sein de $U_2$, la probabilité d'observer un chat est de $1/3$ pour chaque chat. La probabilité d'inclusion de premier ordre du chat échantillonné dans $s_2$ est donc $1/9$.

### 1.1.3 Estimation

Soit $T_y = \sum_{k \in U} y_k$ le total de la variable d'intérêt $y$ sur $U$. Seuls des estimateurs $\hat{T}_y$ du total $T_y$ sont détaillés dans cette section. Ce choix est motivé par deux points : d'une part, les paramètres d'intérêt à La Poste sont des totaux, et d'autre part, les estimateurs de la moyenne $\bar{Y} = T_y/N$ se calculent facilement à partir de $\hat{T}_y$ si $N$ est connue :

$$ \hat{Y} = \hat{T}_y/N. $$

L'objectif d'un sondage est de trouver une bonne stratégie, c'est à dire plan de sondage et estimateur, pour estimer un paramètre d'intérêt, le total dans notre cas. Soit $\hat{T}_y(s)$ la valeur prise par l'estimateur sur l'échantillon $s \in S$. On désigne par $E[\hat{T}_y] = \sum_{s \in S} p(s) \hat{T}_y(s)$ l'espérance sous le plan d'un estimateur $\hat{T}_y$. La qualité d'un estimateur est généralement mesurée par son biais

$$ \text{Biais}(\hat{T}_y) = E[\hat{T}_y] - T_y, $$

sa variance

$$ \text{Var}(\hat{T}_y) = E[(\hat{T}_y - E[\hat{T}_y])^2] $$

et l'**Erreur Quadratique Moyenne (EQM)** qui en découle

$$ E[(\hat{T}_y - T_y)^2] = \text{Var}(\hat{T}_y) + \text{Biais}^2(\hat{T}_y). $$

La précision d'un estimateur est donnée par sa variance. Si la variance de l'estimateur $\hat{T}_y$ est faible sous le plan de sondage $p(\cdot)$ utilisé, alors cette stratégie d'échantillonnage sera performante.

Dans cette section je présente les estimateurs les plus courants, qui sont nécessaires à la compréhension de mes travaux.

#### Estimateur de Horvitz-Thompson

L'estimateur par expansion, ou estimateur de Horvitz-Thompson (HT) du total $T_y$ est l'estimateur le plus simple, introduit par Narain (1951) et Horvitz and Thompson (1952).
<figure><img src="image_11.png" /><figcaption>FIGURE 1.8 – Expansion d'un individu.</figcaption></figure>

Cet estimateur ne nécessite que la connaissance des probabilités d'inclusion d'ordre un, $\pi_k$, $k \in U$.

Soit un individu $k \in U$, avec une probabilité $\pi_k \neq 0$ d'être sélectionné. Cet individu peut être vu comme le représentant de $1/\pi_k$ individus de la population. On appelle poids de sondage, noté $d_k$, la valeur $1/\pi_k$ pour $k \in U$. Dans la figure 1.8, le chat sélectionné a une probabilité d'inclusion de $1/2$, il représente donc 2 individus de la population qui ont des tailles similaires ; lui même et le chat indiqué par la flèche.

L'estimateur de **HT** étend les observations obtenues pour l'échantillon à l'ensemble de la population, comme montré dans la figure 1.9, en multipliant les valeurs $y_k$ observées pour les individus $k \in s$ par leur poids de sondage $d_k$.

<figure><img src="image_12.png" /><figcaption>FIGURE 1.9 – Expansion de l'échantillon.</figcaption></figure>

L'estimateur de **HT** d'un total $T_y = \sum_{k \in U} y_k$ est

$$ \hat{T}_{HT,y} = \sum_{k \in s} d_k y_k. \qquad (1.1) $$
Comme l'échantillon $s$ est l'ensemble des indicatrices d'inclusion $I_k$, $k \in U$, telles que $I_k = 1$, l'estimateur de **HT** peut se réécrire

$$ \hat{T}_{HT,y} = \sum_{k \in U} I_k d_k y_k. $$

Si $\pi_k \neq 0$ pour tout $k \in U$, l'estimateur de **HT** est sans biais. Comme seules les indicatrices d'inclusion sont aléatoires, on a :

$$ E[\hat{T}_{HT,y}] = E[\sum_{k \in U} I_k d_k y_k] = \sum_{k \in U} E[I_k] d_k y_k = \sum_{k \in U} \pi_k d_k y_k = \sum_{k \in U} y_k = T_y. $$

Sa variance est donnée par :

$$ \begin{aligned} \mathrm{Var}(\hat{T}_{HT,y}) &= \mathrm{Var}(\sum_{k \in U} I_k d_k y_k) = \sum_{k \in U} \sum_{k' \in U} \mathrm{Cov}(I_k, I_{k'}) d_k d_{k'} y_k y_{k'} \\ &= \sum_{k \in U} \sum_{k' \in U} \frac{\pi_{kk'} - \pi_k \pi_{k'}}{\pi_k \pi_{k'}} y_k y_{k'}. \end{aligned} \quad (1.2) $$

La variance étant une somme double sur la population, elle ne peut pas être connue à partir de l'échantillon. On ne connaît donc pas la vraie précision d'un estimateur et il faut l'estimer.

Il existe plusieurs estimateurs pour la variance de l'estimateur de **HT**, mais le plus simple consiste à utiliser un estimateur par expansion :

$$ \widehat{\mathrm{Var}}(\hat{T}_{HT,y}) = \sum_{k \in s} \sum_{k' \in s} \frac{1}{\pi_{kk'}} \frac{\pi_{kk'} - \pi_k \pi_{k'}}{\pi_k \pi_{k'}} y_k y_{k'}. \quad (1.3) $$

Cet estimateur de la variance est sans biais si $\pi_{kk'} \neq 0$ pour tout $k, k' \in U$.

Les résultats présentés ci dessus sont valables quelque soit le plan de sondage. Cependant elles se simplifient si un **SASSR** est utilisé, soit directement, soit dans un plan stratifié ou à deux degrés.

## Plan SASSR

Si un **SASSR** de taille $n$ est utilisé pour tirer l'échantillon, on a :

$$ \hat{T}_{HT,y}^{SAS} = \sum_{k \in s} \frac{N}{n} y_k = N \bar{y}, $$

où $\bar{y}$ est la moyenne de $y$ dans l'échantillon. Sa variance est donnée par :

$$ \mathrm{Var}(\hat{T}_{HT,y}^{SAS}) = N^2 \left(1 - \frac{n}{N}\right) \frac{S_y^2}{n}, $$
où $S_y^2 = (N - 1)^{-1} \sum_{k \in U} (y_k - \bar{Y})^2$. Cette variance est estimée sans biais par :

$$ \widehat{\mathrm{Var}}(\hat{T}_{HT,y}^{SAS}) = N^2 \left(1 - \frac{n}{N}\right) \frac{s_y^2}{n}, $$

où $s_y^2 = (n - 1)^{-1} \sum_{k \in s} (y_k - \bar{y})^2$.

## Plan STASSR

Soit une décomposition de $U$ en $H$ strates de taille $N_h$, $h = 1, \dots, H$. On note $T_{yh} = \sum_{k \in U_h} y_k$ le total de $y$ sur la strate $U_h$. Si un STASSR de taille $n$ avec $H$ échantillons $s_h$ de taille $n_h$, $h = 1, \dots, H$, telles que $n = \sum_{h=1}^H n_h$, est utilisé, alors l'estimateur de HT du total $T_y$ est la somme des estimateurs des totaux $T_{yh}$ :

$$ \hat{T}_{HT,y}^{ST} = \sum_{h=1}^{H} \hat{T}_{HT,yh}^{SAS} = \sum_{h=1}^{H} \sum_{k \in s_h} \frac{N_h}{n_h} y_k = \sum_{h=1}^{H} N_h \bar{y}_h, $$

où $\hat{T}_{HT,yh}^{SAS}$ est l'estimateur SASSR du total $T_{yh}$ et $\bar{y}_h$ est la moyenne de $y$ dans l'échantillon $s_h$. Grâce à l'indépendance du tirage entre les strates, sa variance est donnée par :

$$ \mathrm{Var}(\hat{T}_{HT,y}^{ST}) = \sum_{h=1}^{H} \mathrm{Var}(\hat{T}_{HT,yh}^{SAS}) = \sum_{h=1}^{H} N_h^2 \left(1 - \frac{n_h}{N_h}\right) \frac{S_{yh}^2}{n_h}, $$

où $S_{yh}^2 = (N_h - 1)^{-1} \sum_{k \in U_h} (y_k - \bar{Y}_h)^2$ et $\bar{Y}_h$ est la moyenne de $y$ dans la strate $U_h$. Cette variance est estimée sans biais par :

$$ \widehat{\mathrm{Var}}(\hat{T}_{HT,y}^{ST}) = \sum_{h=1}^{H} \widehat{\mathrm{Var}}(\hat{T}_{HT,yh}^{SAS}) = \sum_{h=1}^{H} N_h^2 \left(1 - \frac{n_h}{N_h}\right) \frac{s_{yh}^2}{n_h}, $$

où $s_{yh}^2 = (n_h - 1)^{-1} \sum_{k \in s_h} (y_k - \bar{y}_h)^2$.

Il est possible de décomposer $S_y^2$ comme suit :

$$ S_y^2 = \frac{1}{N-1} \sum_{h=1}^{H} (N_h - 1) S_{yh}^2 + \frac{1}{N-1} \sum_{h=1}^{H} N_h (\bar{Y}_h - \bar{Y})^2 = S_{y,intra}^2 + S_{y,inter}^2, $$

avec $S_{y,intra}^2 = (1/(N-1)) \sum_{h=1}^H (N_h - 1)S_{yh}^2$ mesurant la variabilité au sein des strates (intra-strates) de la variable $y$, $S_{y,inter}^2 = (1/(N-1)) \sum_{h=1}^H N_h (\bar{Y}_h - \bar{Y})^2$ mesurant la variabilité entre les strates (inter-strates) des moyennes $\bar{Y}_h$, $h = 1, \dots, H$ et $(N-1)S_y^2$ mesurant la variabilité
de y sur la population U. La variance de $\hat{T}_{HT,y}^{SAS}$ dépend de $(N-1)S_y^2$ et celle de $\hat{T}_{HT,y}^{ST}$ dépend de $S_{y,intra}^2$.

Sauf cas particuliers, il n'existe pas de propriété générale assurant un gain de précision du STASSR par rapport au SASSR quand un estimateur de HT est utilisé. Toutefois, le plan de sondage stratifié peut permettre un gain de précision par rapport au plan SASSR si la variabilité intra-strates est faible par rapport à la variabilité inter-strates.

## Plan à deux degrés avec SASSR à chaque degré

Soit une décomposition de U en UP $U_i$, $i \in \{1, \dots, N_I\}$. On note $T_{yi}$ le total et $\bar{Y}_i = T_{yi}/N_i$ la moyenne de y dans l'UP $U_i$, $i \in \{1, \dots, N_I\}$. Si un plan à deux degrés avec SASSR de taille $n_I$ au premier degré et SASSR de taille $n_i$, $i \in \{1, \dots, N_I\}$ au second degré est utilisé, alors l'estimateur de HT du total $T_y$ devient :

$$ \hat{T}_{HT,y}^{deg} = \sum_{U_i \in s_I} \sum_{k \in s_i} \frac{N_I}{n_I} \frac{N_i}{n_i} y_k = \sum_{U_i \in s_I} \frac{N_I}{n_I} \hat{T}_{HT,yi} $$

avec $\hat{T}_{HT,yi} = \sum_{k \in s_i} \frac{N_i}{n_i} y_k$ l'estimateur de HT conditionnel au premier degré du total $T_{yi}$, $i \in \{1, \dots, N_I\}$.

La variance de $\hat{T}_{HT,y}^{deg}$ est donnée par :

$$ \mathrm{Var}(\hat{T}_{HT,y}^{deg}) = V_{PSU} + V_{SSU}, $$

avec $V_{PSU} = N_I^2 \left(1 - \frac{n_I}{N_I}\right) \frac{S_{yI}^2}{n_I}$, $S_{yI}^2 = (N_I - 1)^{-1} \sum_{i=1}^{N_I} \left(T_{yi} - \frac{T_y}{N_I}\right)^2$, $V_{SSU} = \frac{N_I}{n_I} \sum_{i=1}^{N_I} N_i^2 \left(1 - \frac{n_i}{N_i}\right) \frac{S_{yi}^2}{n_i}$ et $S_{yi}^2 = (N_i - 1)^{-1} \sum_{k \in U_i} (y_k - \bar{Y}_i)^2$.

$V_{PSU}$ mesure la variance due au premier degré (Primary Sampling Unit) et $V_{SSU}$ mesure la variance due au deuxième degré (Secondary Sampling Unit).

Soit $\bar{y}_i = \sum_{k \in s_i} y_k/n_i$ la moyenne de y dans l'échantillon $s_i$, $i \in \{1, \dots, N_I\}$. Cette variance est estimée sans biais par :

$$ \widehat{\mathrm{Var}}(\hat{T}_{HT,y}^{deg}) = \hat{V}_{PSU} + \hat{V}_{SSU}, $$

avec $\hat{V}_{PSU} = N_I^2 \left(1 - \frac{n_I}{N_I}\right) \frac{s_{yI}^2}{n_I}$, $s_{yI}^2 = (n_I - 1)^{-1} \sum_{i \in s_I} \left(\hat{T}_{HT,yi} - \left(\frac{\sum_{i \in s_I} \hat{T}_{HT,yi}}{N_I}\right)\right)^2$, $\hat{V}_{SSU} = \frac{N_I}{n_I} \sum_{i \in s_I} N_i^2 \left(1 - \frac{n_i}{N_i}\right) \frac{s_{yi}^2}{n_i}$ et $s_{yi}^2 = (n_i - 1)^{-1} \sum_{k \in s_i} (y_k - \bar{y}_i)^2$.
Estimateur de Hájek

L'estimateur de **HT** reste aléatoire même si la variable $y$ est constante. Si toutes les valeurs $y_k, k \in U$ sont égales à une constante $C$, il semble logique que la valeur estimée pour $T_y$ soit égale à $NC$, quel que soit l'échantillon tiré $s$. Or, la valeur de l'estimateur de **HT** du total $T_y$ dans ce cas est égale à :

$$ \hat{T}_{HT,C} = C \sum_{k \in s} d_k, $$

qui n'est généralement pas égale à $NC$ et est aléatoire.

Ce problème peut être corrigé, si la valeur $N$ est connue, en utilisant l'estimateur de Hájek, introduit par Hájek (1971) et donné par :

$$ \hat{T}_{HA,y} = \frac{N}{\sum_{k \in s} d_k} \sum_{k \in s} d_k y_k. $$

Dans le cas où $y_k = C$ pour tout $k \in U$, l'estimateur de Hájek du total est :

$$ \hat{T}_{HA,y} = C \frac{N}{\sum_{k \in s} d_k} \sum_{k \in s} d_k = NC. $$

Soit $\hat{N}_{HT} = \sum_{k \in s} d_k$ l'estimateur de **HT** du total $N = \sum_{k \in U} 1$. L'estimateur de Hájek peut être réécrit comme :

$$ \hat{T}_{HA,y} = \frac{N}{\hat{N}_{HT}} \hat{T}_{HT,y}. \qquad (1.4) $$

L'estimateur de Hájek du total n'est pas sans biais à cause de l'estimateur $\hat{N}_{HT}$ au dénominateur, et sa variance est plus complexe à calculer pour les mêmes raisons. On peut cependant calculer une variance approchée en utilisant des techniques de linéarisation valables pour $n$ et $N$ grands.

L'estimateur de Hájek est un cas particulier d'estimateur par le ratio (voir Särndal et al. (1989), Tillé (2019)).

Propriétés asymptotiques

Les propriétés asymptotiques des estimateurs utilisés en sondage ne peuvent pas s'étudier comme en statistique classique car la population $U$ est finie. Il faut donc définir un contexte spécifique. Considérons le contexte asymptotique décrit par Isaki and Fuller (1982).

On considère une suite infinie de populations finies $U_t$ de taille $N_t$, $t > 1$ telles que pour tout $t$, $U_t$ est inclue dans $U_{t+1}$, $U_t \subset U_{t+1}$, ce qui implique que pour tout $t$, $N_t \le N_{t+1}$. Dans chaque population $U_t$, on observe une variable $Y_t$ ayant pour valeurs $y_{kt}$, $k \in U_t$. Pour tout
$t$, on définit un plan de sondage $p_t$, avec $\pi_{kt}$ et $\pi_{kk't}$ les probabilités d'inclusion de premier et second ordre associées à $p_t$. Un échantillon $s_t$ de taille $n_t$ est tiré dans $U_t$ selon le plan de sondage $p_t$, tel que pour tout $t$, $n_t < n_{t+1}$. Soit $\theta_t$ la valeur du paramètre d'intérêt pour la population $U_t$ et $\hat{\theta}_t$ un estimateur de $\theta_t$ obtenu à partir de $s_t$ dans la population $t$. Quand $t \to \infty$, alors $N_t \to \infty$ et $n_t \to \infty$. Dans la suite, l'indice $t$ sera omis pour alléger les notations.

Dans cette section, nous allons présenter les propriétés asymptotiques de l'estimateur de **HT** de $T_y$ puis de l'estimateur de Hájek de $T_y$.

On considère les hypothèses suivantes (Breidt and Opsomer, 2000) :

H1 On suppose que $\lim_{N \to \infty} \frac{n}{N} = f \in (0, 1)$.

H2 On suppose que les probabilités d'inclusion vérifient les conditions suivantes :

(a) il existe une constante $\lambda$ telle que $\min_{k \in U} \pi_k \ge \lambda > 0$

(b) il existe une constante $c_1$ telle que $\limsup_{N \to \infty} n \max_{k \ne k' \in U} |\pi_{kk'} - \pi_k \pi_{k'}| \le C_1 < \infty$ avec $C_1 > 0$;

H3 On suppose que $\lim_{N \to \infty} T_y/N < \infty$ et qu'il existe une constante positive $C_2$ telle que, pour tout $N$, $\sum_{k \in U} y_k^2/N \le C_2 < \infty$.

**Proposition 1.1.1.** Si les hypothèses précédentes sont vérifiées, alors il existe une constante positive $C$ telle que

$$nE\left(\frac{1}{N}(\hat{T}_{HT,y} - T_y)\right)^2 \le C,$$

pour tout $n$.

*Preuve de la Proposition 1.1.1.* On a

$$\frac{n}{N^2} \mathrm{E}(|\hat{T}_{HT,y} - T_y|^2) = \frac{n}{N^2} \mathrm{Var}(\hat{T}_{HT,y}) = \frac{n}{N^2} \sum_{k \in U} y_k^2 \left(\frac{1}{\pi_k} - 1\right) + \frac{n}{N^2} \sum_{k \in U} \sum_{k' \ne k \in U} y_k y_{k'} \frac{\pi_{kk'} - \pi_k \pi_{k'}}{\pi_k \pi_{k'}}.$$

Selon H2, on peut majorer $(1/\pi_k) - 1$ par $\lambda^{-1}$ pour tout $k \in U$ et selon H3, on peut majorer $(1/N) \sum_{k \in U} y_k^2$ par $C_2$ :
$$
\begin{align*}
& \frac{n}{N^2} \sum_{k \in U} y_k^2 \left( \frac{1}{\pi_k} - 1 \right) + \frac{n}{N^2} \sum_{k \in U} \sum_{k' \neq k \in U} y_k y_{k'} \frac{\pi_{kk'} - \pi_k \pi_{k'}}{\pi_k \pi_{k'}} \\
& \leq \frac{nC_2}{N\lambda} + n \max_{k \neq k' \in U} |\pi_{kk'} - \pi_k \pi_{k'}| \left( \frac{1}{\lambda^2} \left( \frac{1}{N} \sum_{k \in U} |y_k| \right) \right)^2 \\
& \leq \frac{nC_2}{N\lambda} + n \max_{k \neq k' \in U} |\pi_{kk'} - \pi_k \pi_{k'}| \left( \frac{1}{\lambda^2} \sum_{k \in U} y_k^2 \right) \quad \text{selon l'inégalité de Cauchy-Schwarz} \\
& \leq \frac{nC_2}{N\lambda} + n \max_{k \neq k' \in U} |\pi_{kk'} - \pi_k \pi_{k'}| \frac{1}{\lambda^2} C_2 \\
& \leq C \quad \text{par les hypothèses H1 et H2 (b)}
\end{align*}
$$

avec $C > 0$.

□

**Corollaire 1.1.1.** Si les hypothèses précédentes sont vérifiées, alors $\hat{T}_{HT,y}$ est asymptotiquement sans biais pour $T_y$ :

$$
\lim_{N \to \infty} \frac{1}{N} (E(\hat{T}_{HT,y}) - T_y) = 0,
$$

et consistant pour $T_y$ :

$$
\lim_{N \to \infty} P(|\hat{T}_{HT,y} - T_y| > N\epsilon) = 0 \text{ pour tout } \epsilon > 0.
$$

En particulier,

$$
\frac{1}{N}(\hat{T}_{HT,y} - T_y) = O_p\left(\frac{1}{\sqrt{n}}\right).
$$

Preuve du corollaire 1.1.1. Soit $\epsilon > 0$. En utilisant l'inégalité de Markov, on obtient :

$$
P \left( \frac{1}{N} |\hat{T}_{HT,y} - T_y| \geq \varepsilon \right) \leq \frac{\frac{1}{N} E(|\hat{T}_{HT,y} - T_y|)}{\varepsilon},
$$

et l'inégalité de Cauchy-Schwarz

$$
E \left( \frac{1}{N} |\hat{T}_{HT,y} - T_y| \right) \le \sqrt{\frac{1}{N^2} E (|\hat{T}_{HT,y} - T_y|^2)} \le \frac{C}{\sqrt{n}},
$$

par la Proposition 1.1.1. On obtient donc :

$$
\lim_{N \to \infty} P(|\hat{T}_{HT,y} - T_y| > N\epsilon) = 0 \text{ pour tout } \epsilon > 0
$$

et

$$
\lim_{N \to \infty} \frac{1}{N} (E(\hat{T}_{HT,y}) - T_y) = 0.
$$
car $E(\hat{T}_{HT,y} - T_y)/N \le E|\hat{T}_{HT,y} - T_y|/N$.

□

La proposition suivante permet de linéariser l'estimateur de Hájek $\hat{T}_{HA,y} = \sum_{k \in s}^{N} \frac{d_k}{N} \sum_{k \in s} d_k y_k$ du total $T_y$ pour obtenir une variance approximée et un estimateur de variance.

**Proposition 1.1.2.** Si les hypothèses H1-H4 sont vérifiées, alors l'estimateur de Hájek de $T_y$ peut être linéarisé comme suit :

$$ \frac{1}{N} (\hat{T}_{HA,y} - T_y) - \frac{1}{N} \left( \sum_{k \in s} d_k e_k - \sum_{k \in U} e_k \right) = O_p (1/n), $$

où $e_k = y_k - \bar{Y}$, $k \in U$ avec $\sum_{k \in U} e_k = 0$. On obtient en particulier grâce au corollaire 1.1.1 que

$$ \frac{1}{N} (\hat{T}_{HA,y} - T_y) = O_p (1/\sqrt{n}). $$

Preuve de la Proposition 1.1.2. Soit :

$$ D = \frac{1}{N} (\hat{T}_{HA,y} - T_y) - \frac{1}{N} \left( \sum_{k \in s} d_k e_k - \sum_{k \in U} e_k \right) = \frac{1}{N} \hat{T}_{HA,y} - \frac{1}{N} T_y - \frac{\hat{T}_{HT,y}}{N} + \bar{Y} \frac{\hat{N}_{HT}}{N}. $$

On a :

$$ \begin{aligned} D &= \frac{\hat{T}_{HT,y}}{\hat{N}_{HT}} - \frac{\hat{T}_{HT,y}}{N} - \bar{Y} + \bar{Y} \frac{\hat{N}_{HT}}{N} \\ &= \frac{\hat{T}_{HT,y}}{\hat{N}_{HT}} - \bar{Y} + \bar{Y} \frac{\hat{N}_{HT}}{N} - \frac{\hat{T}_{HT,y}}{N} \frac{\hat{N}_{HT}}{N} \\ &= - \left( \frac{\hat{N}_{HT}}{N} - 1 \right) \left( \frac{\hat{T}_{HT,y}}{\hat{N}_{HT}} - \bar{Y} \right) \end{aligned} $$

Selon le corollaire 1.1.1,

$$ \frac{\hat{N}_{HT}}{N} - 1 = \frac{\hat{N}_{HT} - N}{N} = O_p (1/\sqrt{n}). $$

On peut réécrire

$$ \begin{aligned} \frac{\hat{T}_{HT,y}}{\hat{N}_{HT}} - \bar{Y} &= \left( \frac{1}{\hat{N}_{HT}} - \frac{1}{N} \right) \hat{T}_{HT,y} + \frac{1}{N} (\hat{T}_{HT,y} - T_y) \\ &= \frac{N - \hat{N}_{HT}}{N} \frac{N}{\hat{N}_{HT}} \frac{\hat{T}_{HT,y}}{N} + \frac{1}{N} (\hat{T}_{HT,y} - T_y) \end{aligned} $$

avec $\frac{N - \hat{N}_{HT}}{N} = O_p(1/\sqrt{n})$ et $\frac{1}{N}(\hat{T}_{HT,y} - T_y) = O_p(1/\sqrt{n})$ selon le corollaire 1.1.1,
$\frac{N}{\hat{N}_{HT}} = O_p(1)$ selon les hypothèses H1 et H2a et $\frac{\hat{T}_{HT,y}}{N} = O_p(1)$ selon le corollaire 1.1.1 et l'hypothèse H3. On en déduit que

$$ \frac{\hat{T}_{HT,y}}{\hat{N}_{HT}} - \bar{Y} = O_p(1/\sqrt{n}) $$

et donc que $D = O_p(1/n)$.

Comme $\frac{1}{N}(\sum_{k \in s} d_k e_k - \sum_{k \in U} e_k) = O_p(1/\sqrt{n})$ d'après le corollaire 1.1.1

$$ \frac{1}{N} (\hat{T}_{HA,y} - T_y) = D + \frac{1}{N} \left( \sum_{k \in s} d_k e_k - \sum_{k \in U} e_k \right) = O_p \left( 1/\sqrt{n} \right). $$

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

La Proposition 1.1.2 permet d'approximer la variance de l'estimateur de Hájek $\hat{T}_{HA,y}$ de $T_y$ par la variance de $\hat{T}_{HT,E} = \sum_{k \in s} d_k e_k$. On en déduit le corollaire suivant.

**Corollaire 1.1.2.** Si les hypothèses H1, H2, H3 et H4 sont vérifiées, alors la variance de l'estimateur de Hájek $\hat{T}_{HA,y}$ du total $T_y$ peut être approximée par

$$ AVar(\hat{T}_{HA,y}) = Var(\hat{T}_{HT,E}) = \sum_{k \in U} \sum_{k' \in U} \frac{\pi_{kk'} - \pi_k \pi_{k'}}{\pi_k \pi_{k'}} e_k e_{k'} \quad (1.5) $$

où la variable E prend les valeurs $e_k = y_k - \bar{Y}$ sur U.

Cette variance approximée peut être estimée si $\pi_{kk'} \neq 0$ pour tout $k, k' \in U$ par

$$ \widehat{Var}(\hat{T}_{HT,E}) = \sum_{k \in s} \sum_{k' \in s} \frac{\pi_{kk'} - \pi_k \pi_{k'}}{\pi_{kk'} \pi_k \pi_{k'}} \hat{e}_k \hat{e}_{k'} \quad (1.6) $$

avec $\hat{e}_k = y_k - \hat{T}_{HT,y}/N$.

### 1.1.4 Sondage indirect

La base de sondage de la population d'intérêt n'est pas toujours disponible. Dans ce cas, un tirage direct dans la population d'intérêt est impossible et une solution est d'utiliser un sondage indirect.

Un sondage indirect consiste à utiliser une population $U_F$, dite population "Frame" de taille $N_F$ pour laquelle on dispose d'une base de sondage et qui est reliée à la population d'intérêt $U$, de telle façon que tous les individus de $U$ soient reliés au minimum à un individu de $U_F$. Cette condition permet de s'assurer que tous les individus de $U$ peuvent être
<figure><img src="image_14.png" /><figcaption>FIGURE 1.10 – Echantillons tirés via un plan indirect.</figcaption></figure>

échantillonnés avec une probabilité non nulle. Un échantillon $s_F$ est tiré dans $U_F$ et l'échantillon $s$ dans la population d'intérêt est composé de tous les individus reliés à $s_F$. Il n'est pas nécessaire que les individus de $U_F$ soient tous reliés à au moins un individu de $U$, mais cela permet d'obtenir un échantillonnage plus efficace, car un individu de $U_F$ non relié à $U$ n'augmente pas la taille de l'échantillon $s$ s'il est sélectionné.

Par exemple, dans la figure 1.10, il n'existe pas de base de sondage pour la population $U$ des animaux de compagnie en France. On peut tirer indirectement un échantillon d'animaux de compagnie en utilisant la population des français : un échantillon de français est tiré et leurs animaux de compagnie constituent l'échantillon final. Les liens entre les populations sont indiqués en vert. On vérifie bien que tous les chats français sont reliés à au moins un français, mais un français n'est pas nécessairement relié à un chat. Le français tiré dans $U_F$ est lié à deux chats, qui feront partie de l'échantillon de $U$.

Il est aussi possible d'utiliser un plan à deux degrés (voir section 1.1.2) quand la base de sondage de la population d'intérêt n'est pas disponible. Les individus échantillonnés lors d'un sondage à deux degrés sont des **Unités Primaires (UP)** de la population d'intérêt et ceux échantillonnés lors d'un sondage indirect sont des individus reliés à la population d'intérêt.
Le choix de la méthode utilisée dépend donc de l'information disponible.

Avec un plan de sondage indirect, la probabilité de tirer un individu $k \in U$ est complexe à calculer, car elle dépend des liens entre les deux populations. En effet, la probabilité de tirer $k \in U$ est égale à la probabilité de tirer au moins un individu $i$ de $U_F$ lié à $k$. Les poids de sondage calculés par rapport au plan de sondage dans $U$ sont donc inconnus ou très difficilement calculables, contrairement aux poids de sondage calculés par rapport au plan de sondage dans $U_F$, et les estimateurs présentés en section 1.1.3 ne peuvent pas être utilisés tels quels.

Cependant, Deville and Lavallée (2006) ont proposé une méthode, la Méthode Généralisée de Partage des Poids (MGPP), qui contourne cette difficulté. Soit $\theta_{ik}$ un poids de lien entre $i \in U_F$ et $k \in U$, avec $\theta_{ik}$ prenant une valeur positive si $i$ et $k$ sont liés, 0 sinon. La valeur positive de $\theta_{ik}$ est choisie par le statisticien en fonction des informations disponibles. Soit $\tilde{\theta}_{ik} = \theta_{ik}/(\sum_{i \in U_F} \theta_{ik})$ le poids standardisé associé à $\theta_{ik}$, c'est-à-dire tel que $\sum_{i \in U_F} \tilde{\theta}_{ik} = 1$.

La standardisation permet de réécrire le total d'intérêt $T_y = \sum_{k \in U} y_k$ sur $U$ comme un total $T_{\tilde{y}}$ sur $U_F$ de la variable $\tilde{y}_i = \sum_{k \in U} \tilde{\theta}_{ik} y_k$ :

$$T_y = \sum_{k \in U} y_k = \sum_{k \in U} \left( \sum_{i \in U_F} \tilde{\theta}_{ik} \right) y_k = \sum_{i \in U_F} \sum_{k \in U} \tilde{\theta}_{ik} y_k = \sum_{i \in U_F} \tilde{y}_i = T_{\tilde{y}}. \quad (1.7)$$

L'idée est de réécrire le total $T_y$ comme un total d'une nouvelle variable sur la population $U_F$ dans laquelle on effectue le tirage. On peut donc utiliser l'estimateur de **HT** du total $T_{\tilde{y}}$. L'estimateur de **HT** du total $T_{\tilde{y}}$ sur $U_F$ est appelé **estimateur MGPP** du total $T_y$ sur $U$ :

$$\hat{t}_{y1} = \sum_{i \in s_F} \frac{\tilde{y}_i}{\pi_i} = \sum_{i \in s_F} \frac{1}{\pi_i} \left( \sum_{k \in U} \tilde{\theta}_{ik} y_k \right) = \sum_{k \in U} \left( \sum_{i \in s_F} \frac{\tilde{\theta}_{ik}}{\pi_i} \right) y_k.$$

Les seuls individus $k$ de $U$ pour lesquels la somme $\sum_{i \in s_F} \tilde{\theta}_{ik}/\pi_i$ n'est pas nulle sont ceux appartenant à $s$. L'estimateur **MGPP** du total $T_y$ sur $U$ peut donc s'écrire :

$$\hat{t}_{y1} = \sum_{k \in s} \left( \sum_{i \in s_F} \frac{\tilde{\theta}_{ik}}{\pi_i} \right) y_k. \quad (1.8)$$

Supposons que $i \in s_F$ et $k \in U$ soient liés. La Méthode Généralisée de Partage des Poids consiste à utiliser les poids standardisés pour partager le poids de sondage de $i$ entre tous les individus, échantillonnés ou pas, reliés à $k$.

La figure 1.11 permet d'imager le principe du partage des poids et présente un cas très simplifié. Dans cet exemple, on s'intéresse aux chats français sur une zone géographique
<figure><img src="image_15.png" /><figcaption>FIGURE 1.11 – Exemple de poids partagé.</figcaption></figure>

réduite. La population d'intérêt est composée d'un chat et la population Frame est composée de trois personnes, qui sont toutes liées au chat. Il y a donc trois liens, indiqués en bleu. Le statisticien a décidé d'attribuer un poids de liens standardisé de 1/3 à chaque lien.

On tire un individu dans la population Frame avec les probabilités indiquées. Dans ce cas très simple, le chat est indirectement tiré quel que soit la personne échantillonnée et la probabilité d'observer le chat est de 1. Supposons qu'on tire la personne dans l'ellipse. L'individu tiré a un poids de 3. Comme il y a trois personnes permettant d'échantillonner le chat, ce poids doit être partagé entre les trois liens, c'est-à-dire divisé par trois. On obtient donc un poids de 1 pour le chat français. Si la première (resp. troisième) personne en partant du haut était échantillonnée, on obtiendrait un poids de 2/3 (resp. 2) pour le chat. Pour calculer les poids de sondage, il faut donc connaître le nombre de personnes reliées au chat échantillonné, c'est-à-dire 3. Dans le cas où il y a beaucoup de liens entre la population Frame et les individus de U, ce nombre peut être difficile à obtenir, comme dans le cas de La Poste.

L'estimateur **MGPP** a les propriétés d'un estimateur de **HT**, soit une absence de biais pour estimer $T_y$ et une expression connue de la variance :

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

avec $\Delta_{ii'} = (\pi_{ii'} - \pi_i\pi_{i'})/\pi_i\pi_{i'}$ et $\Delta_i = \Delta_{ii}$. Cette variance dépend du choix de la valeur positive des liens pondérés, mais l'estimateur **MGPP** reste sans biais tant que les poids sont standardisés.

La précision d'un estimateur **MGPP** dépend non seulement des poids choisis, mais aussi et surtout des liens entre populations, qui sont fixés. Deville et Lavallée (2006) posent la
question du choix des poids de liens et cherchent des poids “optimaux” qui minimisent la variance de l’estimateur MGPP. Cependant, cette variance dépend aussi des valeurs prises par la variable d’intérêt et on a donc potentiellement autant de poids optimaux que de variables d’intérêt. Par la suite, on appelera poids optimaux, notés $\tilde{\theta}_{ik}^{opt}$, $i \in U_F$, $k \in U$, les poids standardisés qui minimisent la variance de l’estimateur MGPP pour toute variable d’intérêt, c’est-à-dire que pour toute variable d’intérêt $y_k$, $k \in U$, et tout poids standardisé $\tilde{\theta}_{ik}$, $i \in U_F$, $k \in U$ :

$$ \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \sum_{k \in U} \tilde{\theta}_{ik}^{opt} y_k \sum_{k' \in U} \tilde{\theta}_{i'k'}^{opt} y_{k'} \le \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \sum_{k \in U} \tilde{\theta}_{ik} y_k \sum_{k' \in U} \tilde{\theta}_{i'k'} y_{k'} \quad (1.10) $$

## Liens “Tous-pour-Un”

On appelle liens de type **Tous-pour-Un** (TpU) les liens où un individu de $U_F$ n’est lié qu’à un individu de $U$, mais un individu de $U$ peut être lié à plusieurs individus de $U_F$. Par exemple, un français n’est relié qu’à un chat, mais un chat peut être lié à plusieurs français, comme dans la figure 1.11. Si le plan de sondage utilisé pour tirer dans $U_F$ est un **SASSR** ou un plan de Poisson, Deville et Lavallée (2006) montrent qu’il existe des poids optimaux qui minimisent la variance de l’estimateur **MGPP** quelle que soit la variable d’intérêt si les liens sont de type “TpU” et qui sont donnés :

pour le **SASSR**, par $\theta_{ik} = 1$ si $i \in U_F$ et $k \in U$ sont liés, 0 sinon, et

pour le plan de Poisson, par $\theta_{ik} = \frac{1}{(1 - \pi_i)\pi_i}$ si $i \in U_F$ et $k \in U$ sont liés, 0 sinon.

Ces valeurs correspondent aux poids non standardisés. On appelle **MGPP** optimale et estimateur **MGPP** optimal la **MGPP** obtenue en utilisant les poids optimaux et l’estimateur résultant.

## 1.2 La Poste

La Poste est une société anonyme française assurant le transport d'objets postaux (lettres, colis) en France et à l'étranger. Les ancêtres de La Poste sont La Poste à chevaux, créée par Louis XI en 1477, qui effectuait le transport des messages royaux, et l'office des messagers royaux, qui transportait du courrier de particuliers.

Les premières boites aux lettres apparaissent avec la création de la “petite poste” à Paris en 1758-1760. Elle permet la distribution des lettres intra-muros. Avant, les lettres devaient être récupérées au bureau de poste. Le système est par la suite étendu aux grandes villes du
royaume (Bordeaux, Lyon, Nantes, Rouen, Nancy, Strasbourg, Marseille et Lille).

La révolution française transforme le système postal, et entre autres la protection du secret postal, avec le décret sur le secret et l'inviolabilité des lettres en 1791, puis l'ouverture aux services de transport de personnes ou biens aux particuliers en 1794. Le monopole postal est déclaré en 1801 et la direction générale des Postes est créée en 1804.

Au 18e siècle, La Poste se modernise et s'étend aux campagnes. Le retour à l'envoyeur, l'utilisation de bateaux à hélice ou de bureaux de poste ambulants, permettant le tri à bord de wagons des chemins de fer, les timbres-poste, l'uniforme et la carte postale sont créés sur cette période. Les femmes sont engagées à La Poste depuis 1877.

Au 20e siècle, l'essor continue avec la création et l'âge d'or de l'Aéropostale, la création de La Poste automobile rurale (qui remplace les bicyclettes), les distributeurs automatiques de timbres, les premières machines à affranchir... et les premières grèves.

Le 21e siècle voit l'ouverture progressive du marché postal à la concurrence. La Poste change de statut pour devenir SA en 2010.

L'histoire de La Poste est présentée en détail au Musée de La Poste (34 Bd de Vaugirard, 75015 Paris) ou sur leur site https://www.museedelaposte.fr/fr/histoire-de-la-poste.

## 1.2.1 Plans de sondage SYCI 1 et 2

La Poste distribue plusieurs types de produits (lettres recommandées, lettres vertes, co-lis...), appelés objets postaux, ayant chacun des caractéristiques différentes, comme les dimensions de l'enveloppe, le poids, le type et le prix du timbre ou encore le délai maximum de livraison. Les objets postaux comprennent les communications (lettres...) et les objets de marchandise (colis...) et sont livrés lors des tournées de facteur. Une tournée est un ensemble d'adresses postales pour lesquelles un facteur va livrer du courrier.

La régulation de La Poste dans le cadre du Service Universel (SU) est assurée par l'AR-CEP (Autorité de Régulation des Communications Électroniques, des Postes et de la distribution de la presse) qui requiert la mise à disposition de chiffres précis dans le domaine postal afin d'effectuer une régulation efficace. La Poste cherche donc à avoir une vision globale de la mise en distribution des objets. Entre autres, elle souhaite connaître le volume et la composition du trafic postal mensuel, trimestriel et annuel. La Poste s'intéresse aussi au traitement des objets et à l'organisation des tournées, en regardant notamment les délais de livraison, si un objet a été trié par machine de tri et le mode de transport utilisé pour la livraison.

La Poste ne pouvant pas observer toutes les tournées distribuées en France Métropolitaine, elle a mis en place depuis 1994 une étude par sondage, SYCI 1, (SYstème de Collecte
de l'Information), afin d'estimer la composition du trafic et assurer un suivi du traitement des objets postaux. SYCI est un observatoire des flux d'objets distribués et de l'organisation de la distribution.

Dans cette section, nous allons présenter le plan de sondage SYCI 1 utilisé entre 1994 et 2015. Des réorganisations des centres de tri ont poussé La Poste a revoir sa méthodologie en 2015. Nous présenterons le plan initialement prévu en 2015 et enfin le plan actuel SYCI 2.

### SYCI 1

La population d'intérêt est composée des objets distribués un trimestre donné. Par la suite, on simplifiera cette population en se restreignant à un mois donné. La Poste cherchant aussi à observer l'organisation de la distribution, elle a décidé d'observer les tournées, aussi appelées SOrties-Jours (SOJ) ou tournées-jours.

Les tournées peuvent être vues comme un groupe d'objets, de taille trop importante pour que tous les objets d'une tournée puissent être observés. L'étude historique SYCI 1 consistait donc à tirer avec un plan à deux degrés (voir section 1.1.2) un échantillon de tournées puis un échantillon d'objets au sein de chaque tournée échantillonnée. Les objets observés sont appelés objets-échantillon.

Les plans de sondage utilisés pour tirer les tournées puis les objets au sein des tournées sont des STASSR (voir section 1.1.2). Le poids de sondage final d'un objet est donc le poids de sondage de sa tournée d'appartenance multiplié par son poids de sondage au sein de la tournée.

Ce plan de sondage repose sur l'hypothèse que la population des tournées est suffisamment stable pour que la base de sondage corresponde à la réalité du terrain. Cependant, les centres de distribution se réorganisent régulièrement pour permettre d'adapter la distribution selon le volume de courrier. Ces réorganisations ont rendu inadéquate la base de sondage des tournées. La Poste a donc mis en place en 2015 un nouveau plan de sondage, SYCI 2, qui modifie le tirage des tournées, mais garde le plan à deux degrés.

### Plan SYCI 2 initialement prévu

Pour tirer les tournées, plusieurs solutions existent. On peut par exemple effectuer un tirage à deux degrés en tirant des bureaux puis des tournées au sein des bureaux. Cependant, cette méthode nécessite un nombre important d'enquêteurs et de moyens.

La Poste a décidé d'échantillonner les tournées via un sondage indirect (voir section 1.1.4), en utilisant la population des adresses-jours, qui est accessible, comme population Frame et la population des tournées comme population cible. Ce choix est motivé par la stabilité de
la population des adresses, qui n'est pas sensible aux réorganisations des bureaux. On tire un échantillon d'adresses et un échantillon de jours, puis les adresses-jours sont formées en associant aléatoirement un jour à une adresse. Le jour donné, l'enquêteur va voir le bureau qui s'occupe de l'adresse et repère la tournée qui dessert l'adresse tirée. Par simplicité, la population des adresses-jours sera désignée comme population des adresses dans ce qui suit.

<figure><img src="image_16.png" /><figcaption>FIGURE 1.12 – MGPP initialement prévue par La Poste.</figcaption></figure>

L'utilisation d'un sondage indirect va de pair avec l'utilisation d'une MGPP. Or, pour pouvoir utiliser la MGPP, il faut pouvoir calculer les poids standardisés (voir section 1.1.4). Pour cela il est nécessaire d'observer tous les liens entre la tournée échantillonnée et la population Frame, ce qui, pour La Poste, revient à pouvoir lister toutes les adresses de la tournée. Une tournée étant composée de 500 adresses en moyenne, comme indiqué dans la figure 1.12, récupérer toutes les adresses prend du temps. De plus, les adresses ne sont pas toujours correctement indiquées, ce qui rend leur identification complexe et nécessite l'aide du facteur. Le facteur ne pouvant pas prendre de retard pour répondre aux questions de l'enquêteur, il n'est généralement pas possible de récupérer les adresses d'une tournée. L'utilisation du sondage indirect, et de la MGPP associée, pose donc problème.

## Plan SYCI 2

La solution adoptée par La Poste consiste à utiliser une deuxième population intermédiaire, liée à la fois à la population des adresses et des tournées. La population en question est la population des cases de tri, qui sont les cases permettant de trier les objets d'une tournée. Cette population est naturellement entre la population des adresses et la population des tournées.

Pour comprendre les choix de La Poste, regardons le fonctionnement d'un bureau de Poste avant le départ des facteurs en tournée. Tous les matins (en jour ouvré), les objets
<figure><img src="image_17.png" /><figcaption>FIGURE 1.13 – Un casier (à gauche) et une case de tri (à droite).</figcaption></figure>

postaux à distribuer sont livrés dans les centres de tri. Ces objets sont ensuite triés dans des cases, appelées cases de tri, qui sont attribuées à une ou plusieurs adresses. Les objets à destination de ces adresses sont classés dans la case correspondante. Ces cases peuvent être de différents types, comme indiqué dans la figure 1.14, mais le principe reste le même quel que soit le type. Les cases sont regroupées en casiers similaires à celui de la figure 1.13. Les

Case :

CHM

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

CM

<figure><img src="image_19.png" /><figcaption>FIGURE 1.14 – Différents types de cases.</figcaption></figure>

tournées sont formées à partir des casiers. Pour simplifier, on supposera par la suite qu'une tournée est composée de tous les objets contenus dans un casier et qu'un casier contient seulement les objets de la tournée correspondante.

Pour tirer un échantillon de tournées via la population des cases, deux tirages indirects (voir section 1.1.4) sont réalisés. Le premier permet d'obtenir un échantillon de cases à partir de la population des adresses et le second permet d'obtenir un échantillon de tournées à partir des cases (indirectement) tirées. On parle de sondage indirect double. L'échantillon final est composé des tournées liées aux cases indirectement tirées, qui sont elles-mêmes liées aux
<figure><img src="image_20.png" /><figcaption>FIGURE 1.15 – MGPP double utilisée à La Poste.</figcaption></figure>

adresses sélectionnées, comme indiqué dans la figure 1.15.

La population intermédiaire n'influe pas sur les liens entre la population Frame et la population cible. Une adresse fait toujours partie de la même tournée, que sa case de tri soit observée ou non. Tirer un échantillon par sondage indirect double revient donc à tirer un échantillon par sondage indirect, où les liens sont observés en deux temps plutôt que directement. Par la suite, le sondage indirect avec des liens observés en un seul temps sera appelé sondage indirect simple. Pour déterminer la tournée d'appartenance d'une adresse, il faut identifier sa case de tri puis la tournée, plutôt qu'identifier directement la tournée. Cette complexité d'observation des liens permet d'obtenir des estimateurs dont l'expression est plus complexe que ceux obtenus avec un sondage indirect simple.

L'estimateur est obtenu en appliquant une MGPP pour chaque ensemble de liens (voir section 1.1.4). On appelle donc cette méthode une MGPP double. Un premier ensemble de poids de liens est choisi pour les liens entre les adresses et les cases, et un second entre les cases et les tournées. Pour standardiser les poids de liens, il faut donc identifier toutes les adresses des cases échantillonnées et toutes les cases des tournées échantillonnées, comme indiqué dans la figure 1.15. Les liens rouges sont échantillonnés et pour pouvoir standardiser leur poids, il faut observer les liens verts, ce qui revient à observer 10 adresses pour une case et 50 cases pour une tournée en moyenne. Une tournée livrant 500 adresses en moyenne, cette méthode permet de réduire à 60 le nombre de liens observés par rapport à la MGPP initialement prévue, et est donc plus facile à mettre en place.

La MGPP double est un cas particulier de MGPP, comme on le verra dans la section 1.3. Cependant, la MGPP double utilisée à La Poste cause une perte de précision importante comparée aux estimateurs directs obtenus avec SYCI 1, comme le montre la figure 1.16. Sur cette figure, les estimations (normalisées pour des raisons de confidentialité) du trafic total sont indiquées en noir et les intervalles de confiance des estimateurs sont en rouge. On voit
Indicateur Normalisé "Nombre d'objets Courrier" - de Janvier 2009 à février 2023

Estimation SYCI - avec bandes de confiance

<figure><img src="image_21.png" /><figcaption>FIGURE 1.16 – Perte de précision des estimateurs de La Poste.</figcaption></figure>

qu'en 2015, année de transition de **SYCI** 1 à **SYCI** 2, les intervalles de confiance s'élargissent.
La **MGPP** n'ayant pas pu être mise en place, la comparaison entre les estimateurs **MGPP** et **MGPP** double ne peut pas être réalisée dans le cas de La Poste.

Les travaux menés lors de ma thèse s'inscrivent dans ce contexte de perte de précision
due à la transition de **SYCI** 1 à **SYCI** 2. De plus, La Poste prévoit de réduire les tailles
d'échantillons dans un futur proche, ce qui va dégrader la précision des estimateurs. Une
partie de mon travail porte sur l'amélioration de la **MGPP** double pour améliorer la précision
des estimateurs et limiter l'impact de la réduction de la taille des échantillons.

Dans mes travaux sur la **MGPP** (chapitre 2 et chapitre 3), le tirage à deux degrés des
objets est omis. On considère par simplicité que la population d'intérêt est la population des
tournées et que la variable d'intérêt est le nombre d'objets postaux distribués par une tournée.
Le tirage s'effectue par sondage indirect double en utilisant la population des adresses comme
population Frame.

La Poste effectue un travail important de traitement des données récoltées pour s'assurer
de leur fiabilité. De plus, La Poste utilise des méthodes de repondération des estimateurs,
comme le calage (Deville and Särndal, 1992) ou la winsorisation (Kokic and Bell, 1994).
Certains de ces traitements, notamment la winsorisation, biaisent les estimateurs. Bien que le biais causé par ces traitements soit une source de préoccupation pour La Poste, ce problème n'entre pas dans le cadre de la thèse et on étudiera dans la suite les estimateurs avant leur traitement par La Poste.

## 1.2.2 Données auxiliaires

La Poste possède des bases de données de diverses provenances, contenant des informations sur les populations utilisées pour le tirage : les adresses, les cases, les tournées et les objets. Dans cette section, on présente trois sources de données différentes : les données récoltées via *SYCI 2*, qui sont actuellement utilisées pour les estimations, les données permettant d'obtenir la base de sondage, dont une grande partie est inutilisée à l'heure actuelle, et les données récoltées par les machines de tri, dont l'utilisation fait partie de mes sujets de recherche et est actuellement testée à La Poste.

### Données récupérées dans *SYCI 2*

Les données présentées ici sont les données récoltées par les enquêteurs dans le cadre de l'étude *SYCI 2*. Les données récoltées forment trois tables séparées : une table d'information sur les objets observés, une table sur les tournées observées et une table sur les cases observées.

La figure 1.17 montre une fiche de saisie objet. Les informations récoltées par les enquêteurs comprennent, entre autres, le type de produit (chronopost, lettre verte...), la date de dépôt de l'objet, le montant d'affranchissement, le poids, les dimensions et l'épaisseur du pli, la nature du contenant (film plastique pour les pubs ou revues, enveloppe à fenêtre...) et le contenu (s'il y a un objet, du papier...). Ces données composent la table des objets observés. En pratique, elles sont utilisées pour calculer les variables d'intérêt, qui sont souvent des regroupements de types de produit par exemple. Pour simplifier, on supposera par la suite que les variables d'intérêt sont directement récoltées via le formulaire 1.17.

La figure 1.18 montre une fiche de saisie *SOJ*. Il s'agit des informations récoltées sur la tournée, qui permettent de s'assurer du bon déroulement de l'enquête. La strate d'apparition de la tournée et la date d'observation sont indiquées, ainsi que des indications sur l'état de la tournée qui peuvent affecter la saisie des informations. Par exemple, si la tournée n'a pas pu être observée entièrement.

Des informations permettant d'effectuer un suivi de l'organisation de la distribution des objets postaux sont aussi récoltées. Par exemple, le préparateur de la tournée n'est pas toujours celui qui la distribue et le site de préparation de la tournée n'est pas toujours le site de distribution. Ces informations sont donc récoltées, en plus de la période et du mode de
<figure><img src="image_22.png" /><figcaption>FIGURE 1.17 – Fiche de description des objets.</figcaption></figure>

livraison.

Afin de pouvoir mettre en place une **MGPP** double, La Poste doit observer les liens entre les populations. Cette information est récupérée via la fiche 1.19. L'enquêteur indique les cases composant le casier de la tournée. Cette fiche permet aussi de récupérer des informations au niveau des cases, puisque, pour chaque case du casier, le type de case et ses dimensions sont récoltées.

La table des données sur les tournées contient les informations récoltées via la fiche 1.18. On y ajoute le nombre de cases dans la tournée, obtenu grâce à la fiche 1.19. Les autres informations collectées via 1.19 composent la table des données sur les cases.

Ces données sont vérifiées et validées après récolte. De plus, les enquêteurs suivent une formation afin de s'assurer de la fiabilité de leur captation. Les données récoltées via **SYCI 2** sont donc considérées comme fiables.
<figure><img src="image_23.png" /><figcaption>FIGURE 1.18 – Fiche de description des tournées.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>FIGURE 1.19 – Fiche de description des liens entre cases et tournées.</figcaption></figure>

Référentiel des Adresses Organisées (RAO)

La table RAO contient des informations sur les adresses occupées en France Métropoli-
taine. Cette table, croisée avec les données de l'INSEE, permet d'obtenir une base de sondage
pour la population des adresses. Cependant, elle contient des informations sur les tournées
et cases d'appartenance des adresses, comme on peut le voir dans la figure 1.20.

Pour une adresse, RAO connaît son type (rue, numéro de rue, commune...) mais donne
aussi l'identifiant de la tournée qui dessert cette adresse ainsi que l'identifiant de la case
contenant les objets à destination de cette adresse. On remarque aussi que RAO permet de
rechercher des informations directement sur un casier ou une tournée.

Pour chaque tournée et case, RAO contient la liste des adresses la composant. Malheu-
reusement, ces informations ne prennent pas en compte les réorganisations journalières et
ponctuelles des centres de tri et sont parfois manquantes. Les tournées et cases indiquées
par RAO pour une adresse, si elles existent, ne sont généralement pas les tournées et cases
observées sur le terrain.

Cependant, RAO reste une source importante de données sur les cases. Dans ma thèse,
je propose une façon d'utiliser ces données pour améliorer la MGPP double actuellement
utilisée à La Poste.
<figure><img src="image_25.png" /><figcaption>FIGURE 1.20 – Aperçu des données de RAO.</figcaption></figure>

## Traitement Automatisé de l'Enveloppe (TAE)

TAE est la table composée des données récoltées par machine de tri. Les machines prennent une photo de l'objet trié, puis utilisent des algorithmes de reconnaissance d'images pour lire les inscriptions de l'enveloppe.

A partir des photos, les machines récoltent plusieurs informations : l'adresse de destina-
tion, le type d'objet, la date et le montant d'affranchissement, la machine à affranchir utilisée
(si utilisée) et la zone de départ. Ces informations sont récoltées dans le respect du secret
des correspondances et du Règlement Général sur la Protection des Données. Les machines
enregistrent aussi des informations sur le traitement, comme les dates de passage en machine
de tri, la photo ou encore l'identifiant imprimé sur le contenant de l'objet (souvent une en-
veloppe) par la machine de tri. Cet identifiant est imprimé sur le contenant par les machines
sous forme de code barre, similaire à celui montré dans la photo 1.22.

La description des objets réalisée par les machines est basée sur une photo Noir et Blanc de
l'objet, comme on peut le voir dans la figure 1.21. De plus, la machine de tri ne peut prendre
qu'une photo : recto ou verso. Les informations décrites dans le paragraphe précédent sont
inscrites sur le recto et ne peuvent donc pas être obtenues par les machines pour tous les
objets. Par exemple, une photo du verso d'un objet ne contient pas d'information sur le
type, l'adresse de destination, l'affranchissement... Mais elle peut contenir des informations
sur l'expéditeur. Les timbres "marianne", qui se reconnaissent à leur couleur (vert, rouge ou
gris) ne sont pas facilement différentiables si la photo est en noir et blanc, comme montré
<figure><img src="image_26.png" /><figcaption>FIGURE 1.21 – Photos prises par les machines de tri.</figcaption></figure>

dans la figure 1.23. Le prix du timbre, voire le type de pli, n'est alors pas reconnu par la machine. Ce dernier problème est voué à disparaître avec la suppression des timbres marianne rouge et noir.

D'autres informations, récupérées dans **SYCI 2**, ne peuvent pas être lues par les machines. Entre autres, le poids, l'épaisseur, le contenu... Ce qui s'explique par la différence d'objectif de **TAE** et **SYCI**. Les informations collectées dans **TAE** permettent de distribuer l'objet à la bonne adresse et dans les délais impartis, alors que les données collectées dans **SYCI** permettent d'estimer le trafic postal.

**TAE** est moins fiable que **SYCI 2**, de part sa sensibilité à la prise de photo. Une photo floue, trop sombre, trop claire, avec des objets agglomérés rendent la reconnaissance des informations par les machines erronée. A cela s'ajoutent les problèmes liés aux traitements et au transfert des données, qui peuvent être mal transmises d'une machine à une autre, voire perdues lors des coupures réseaux. Dans ce dernier cas, on observe dans **SYCI 2** un objet, traité via machine de tri, mais absent de **TAE**. **TAE** souffre aussi de problème de doublons dû à la gestion des données, car un identifiant objet peut être réutilisé plusieurs fois.

**TAE** et **SYCI** récupèrent tous les deux des variables similaires, mais il est impossible de remplacer **SYCI** par **TAE** à l'heure actuelle pour observer le trafic, à cause des problèmes mentionnés dans cette section (données incomplètes et erronées). L'objectif principal de ma thèse est de trouver un moyen d'utiliser **TAE**, malgré ses inconvénients, pour améliorer la
<figure><img src="image_27.png" /><figcaption>FIGURE 1.22 – Code barre imprimé par la machine de tri.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIGURE 1.23 – Timbres Marianne en couleur (haut) et en noir et blanc (bas).</figcaption></figure>

précision des estimations réalisées via SYCI 2.

## 1.3 Contributions

L'objet de ma thèse est d'améliorer la précision des estimateurs du trafic postal. Par souci de simplicité nous ne prendrons pas en compte les procédures de calage ni de winsorisation utilisées par La Poste (voir section 1.2.1) dans l'étude des estimateurs.

L'étude de la précision des estimateurs considérés s'effectue selon deux axes. D'un côté, nous cherchons à réduire la variance des estimateurs actuels SYCI 2, considérés sans biais, notamment via des travaux sur la MGPP (voir section 1.1.4) et la MGPP double (voir section 1.2.1). D'un autre côté, nous développons des méthodes et des estimateurs qui utilisent
comme information auxiliaire les bases de données massives disponibles à La Poste. Dans cette section, je présente rapidement les travaux effectués pendant ma thèse. Dans l'ensemble de mes travaux de thèse, on ne considère que l'estimateur du trafic postal mensuel, qui est le total des objets postaux distribués en France pendant un mois donné.

1.3.1 Sondage indirect et MGPP

Cette section porte sur la réduction de la variance de l'estimateur MGPP du trafic total mensuel actuellement utilisé à La Poste. Je présente très brièvement les résultats principaux publiés dans l'article "*Many-to-One indirect sampling with application to the French postal traffic estimation*" (Annals of Applied Statistics (AoAS), 2023), dont le contenu figure dans le chapitre 2, ainsi que les résultats présentés dans le chapitre 3, "*Optimal Weights for double Many-To-One Generalized Weight Share Method*".

Chaptitre 2 : *Many-to-One indirect sampling with application to the French
postal traffic estimation (AoAs)*

On a vu en section 1.2.1 que l'utilisation d'une MGPP double a diminué la précision des
estimateurs SYCI 2, utilisés actuellement par La Poste, comparé à l'estimateur de HT SYCI
1 calculé à partir d'un échantillon direct de tournées et utilisé jusqu'en 2015 par La Poste.
Dans cet article, je cherche à comprendre l'origine de cette perte de précision.

Pour faciliter la lecture, la MGPP définie en section 1.1.4 sera appelée MGPP simple et on dira qu’un estimateur du total $T_y$ est direct s’il est obtenu à partir d’un échantillon sélectionné directement dans la population cible $U$. Dans l’article, la population d’intérêt est notée $U_T$ (target population) mais la notation $U$ sera conservée dans ce chapitre par souci de simplicité.

Dans un premier temps, je me suis penchée sur les conditions d'existence des poids de
liens optimaux définis par l'équation (1.10), c'est à dire les poids de liens qui permettent de
minimiser la variance, donnée par l'équation (1.9), de l'estimateur MGPP simple du total
$T_y$ pour une variable d'intérêt y quelconque. Comme mentionné en section 1.1.4, Deville and
Lavallée (2006) cherchent des poids optimaux pour les plans SASSR et de Poisson, pour
tout type de liens. Ils ont montré que ceux-ci existent si et seulement si les liens entre les
populations sont de type "Tous-pour-Un" (TpU). Dans l'article Medous et al. (2023a), je
montre que l'existence de poids optimaux pour la MGPP simple est une propriété des liens
de type TpU qui peut être étendue à tous les plans de sondage satisfaisant la Δ-propriété
définie ci-dessous.

Soit $N_F$ la taille de la population Frame $U_F$ et, pour $k \in U$, $U_{Fk}$ la sous-population de
$U_F$ de taille $N_{Fk}$ composée des individu reliés à $k$. On note

$$ \Delta_{ii'} = \frac{\pi_{ii'} - \pi_i \pi_{i'}}{\pi_i \pi_{i'}}, \quad i, i' \in U_F, $$

$\Delta = (\Delta_{ii'})_{i,i' \in U_F}$ la matrice de taille $N_F \times N_F$ et $\Delta_{kk'} = (\Delta_{ii'})_{i \in U_{Fk}, i' \in U_{Fk'}}$ la sous-matrice de $\Delta$ de taille $N_{Fk} \times N_{Fk'}$ composée des éléments en position $i, i'$ si $i$ (resp. $i'$) est lié à $k$ (resp $k'$) et $\Delta_k = \Delta_{kk}$. Dans la figure 1.24, l'individu $k$ de la population cible $U$ est lié aux individus notés 1, 2 et 3 dans la population $U_F$. La sous-matrice $\Delta_k$ est donc composée des 9 valeurs $\Delta_{ii'}$, $i, i' = 1...3$, situées en haut à gauche de la matrice (en rouge). On peut similairement obtenir les trois autres matrices de la figure 1.24.

<figure><img src="image_29.png" /><figcaption>FIGURE 1.24 – Décomposition de la matrice $\Delta$ dans un cas TpU.</figcaption></figure>

Quand les liens sont de type TpU, les sous-populations $U_{Fk}$ forment une partition de $U_F$ et on suppose les individus de $U_F$ ordonnés selon $U_{Fk}, k \in U$, (voir le graphique de gauche de la figure 1.24). On a alors (voir le graphique de droite de la figure 1.24)

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

Soit $1_k$ le vecteur de taille $N_{Fk}$ dont les composantes sont égales à 1. Si les liens sont de type TpU, alors un plan de sondage satisfait la $\Delta$-propriété si

— pour tout $k \in U$, $\Delta_k$ est inversible,

— pour $k \neq k' \in U$,

$$ \Delta_{kk'} = c_{kk'} 1_k 1_{k'}^t \quad \text{avec } c_{kk'} \text{ ne dépendant ni de } i \text{ ni de } i'. \quad (1.11) $$

Les plans de Poisson, *Sondage Aléatoire Simple Sans Remise (SASSR)* et, sous certaines conditions, *SASSR stratifié (STASSR)* satisfont la $\Delta$-propriété (voir Medous et al., 2023a).
Je montre dans Medous et al. (2023a) que, si les liens sont de type TpU et que le plan de sondage dans $U_F$ satisfait la $\Delta$-propriété, il existe alors des poids de liens standardisés $\tilde{\theta}_{ik}^{opt}$, $i \in U_F$, $k \in U$, donnés par :

$$ (\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{pour tout } k \in U_T, \quad (1.12) $$

qui minimisent la variance de l'estimateur MGPP simple quelle que soit la variable d'intérêt.

La $\Delta$-propriété permet de simplifier l'expression de la variance de l'estimateur par MGPP simple et on montre dans ce cas que la différence entre la variance de l'estimateur MGPP simple $\hat{t}_{y1}$ utilisant des poids quelconques $\tilde{\theta}_{ik}$, $i \in U_F$, $k \in U$, et la variance de l'estimateur MGPP simple optimal $\hat{t}_{y1}^{opt}$ utilisant les poids optimaux $\tilde{\theta}_{ik}^{opt}$ est donnée par :

$$ \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}}), \quad (1.13) $$

avec $\hat{t}_{\tilde{\theta}_k^{opt}} = \sum_{i \in s_F} d_i \tilde{\theta}_{ik}^{opt}$ l'estimateur de HT du total $t_{\tilde{\theta}_k^{opt}} = \sum_{i \in U_F} \tilde{\theta}_{ik}^{opt} = \sum_{i \in s_F} d_i (\tilde{\theta}_{ik} - \tilde{\theta}_{ik}^{opt})$ l'estimateur de HT du total $t_{\tilde{\theta}_k - \tilde{\theta}_k^{opt}} = 0$, pour tout $k \in U_T$. L'équation (1.13) montre que la différence entre la variance de l'estimateur MGPP simple quelconque et la variance de l'estimateur MGPP simple optimal du total $T_y$ dépend de la variance de la différence entre les poids standardisés $\tilde{\theta}_k - \tilde{\theta}_k^{opt}$, $k \in U$. La standardisation dépend de la répartition des liens entre les populations et influence aussi l'augmentation de variance liée à l'utilisation de poids quelconques, comme illustré dans les études par simulation (voir section 2.4).

Comme discuté dans la section 1.2.1, le premier plan indirect prévu par La Poste prévoyait l'utilisation d'une MGPP simple optimale. La Poste souhaitait connaître la performance, en terme de variance, de l'estimateur MGPP simple optimal, comparé à l'estimateur de HT direct de $T_y$. Il n'est en général pas possible de déterminer lequel des deux estimateurs sera le plus précis en terme de variance (Kiesl, 2016). Dans l'article, je me penche sur un cas particulier où un plan de Poisson est utilisé pour tirer l'échantillon direct $s$ dans $U$ et l'échantillon $s_F$ dans $U_F$, avec $\pi_i$ (resp. $\pi_k$) la probabilité d'inclusion d'ordre un de l'individu $i \in U_F$ (resp. $k \in U$) utilisée pour tirer l'échantillon $s_F$ (resp. l'échantillon direct $s$). La probabilité $\pi_k$, $k \in U$ est obtenue comme suit :

$$ \pi_k = P(\text{au moins un } i \in U_F \text{ lié à } k \text{ est échantillonné dans } s_F) = 1 - \prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}. $$

Dans ce cas, je montre que la variance de l'estimateur de HT direct de $T_y$ est plus faible que la variance de l'estimateur MGPP simple optimal de $T_y$ (voir Proposition 2.2.3).

Pour La Poste, $U_F$ est la population des adresses et $U$ celle des tournées. Ces deux
populations sont reliées par des liens de type TpU et La Poste utilise un STASSR tel que toutes les adresses desservies par une même tournée appartiennent à la même strate. Le plan de sondage utilisé à La Poste satisfait la Δ-propriété et il existe donc, en théorie, des poids de liens qui permettent de minimiser la variance de l'estimateur MGPP simple initialement prévu par La Poste (voir section 1.2.1) quelle que soit la variable d'intérêt. Ces poids sont égaux à l'inverse du nombre d'adresses par tournées, qui vaut 500 en moyenne. Comme mentionné en section 1.2.1, l'utilisation de ces poids optimaux ne peut pas être effectuée car il n'est pas possible d'identifier, ni de compter, toutes les adresses d'une tournée.

La Poste a donc été contrainte de réduire le nombre de liens à observer. Pour ce faire, elle a utilisé la MGPP double décrite dans la section 1.2.1. Mon second objectif dans l'article Medous et al. (2023a) est d'introduire le principe général de la MGPP double et d'étudier ses propriétés.

Soit une population intermédiaire $U_M$, reliée à $U_F$ et à $U$ de telle manière que tout individu $k \in U$ soit relié à au moins un individu de $U_M$ et que tout individu de $U_M$ soit relié à au moins un individu de $U_F$. Dans le cas de La Poste, $U_M$ est la population des cases de tri. On attribue à chaque lien entre $i \in U_F$ et $j \in U_M$ (resp. entre $j \in U_M$ et $k \in U$) un poids de liens $\theta_{ij}^{FM}$ (resp. $\theta_{jk}^{MT}$) tel que $\theta_{ij}^{FM} > 0$ si $i$ et $j$ sont liés (resp. $\theta_{jk}^{MT} > 0$ si $j$ et $k$ sont liés), 0 sinon. Le lien entre $i$ et $k$ est pondéré par $\theta_{ik} = \sum_{j \in U_M} \theta_{ij}^{FM} \theta_{jk}^{MT}$ et la MGPP double est donc un cas particulier de MGPP simple. Soit $\tilde{\theta}_{ij}^{FM}$, $\tilde{\theta}_{jk}^{MT}$ et $\tilde{\theta}_{ik}$ les poids standardisés associés à $\theta_{ij}^{FM}$, $\theta_{jk}^{MT}$ et $\theta_{ik} = \sum_{j \in U_M} \theta_{ij}^{FM} \theta_{jk}^{MT}$. Ces poids sont calculés de manière à avoir

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

On parle alors de standardisation globale. L'estimateur MGPP double est obtenu en remplaçant $\tilde{\theta}_{ik}$ par $\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}$ dans l'expression de l'estimateur MGPP simple donnée par l'équation (1.8) :

$$ \hat{t}_{y2} = \sum_{k \in s} \left( \sum_{i \in s_F} \frac{1}{\pi_i} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} \right) y_k. \quad (1.14) $$

On considère le cas où les liens entre les trois populations sont de type TpU-TpU, c'est à dire TpU entre $U_F$ et $U_M$ puis entre $U_M$ et $U$. C'est le cas de La Poste, comme illustré par la figure 1.25. La somme $\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}$ ne contient dans ce cas qu'un seul élément non nul pour un certain $i \in U_F$ et $k \in U$.

Dans l'exemple de la figure 1.25, cet élément non nul est le produit des poids associés aux deux liens rouges. Les poids $\tilde{\theta}_{ij}^{FM}$ et $\tilde{\theta}_{jk}^{MT}$ utilisés à La Poste sont respectivement l'inverse du nombre d'adresses par case et l'inverse du nombre de cases par tournée. Le poids utilisé par La Poste dans l'exemple de la figure 1.25 est donc $1/5 * 1/3 = 1/15$. Dans le cas de La Poste,
<figure><img src="image_30.png" /><figcaption>FIGURE 1.25 – TpU-TpU observé par La Poste.</figcaption></figure>

ce poids est généralement différent du poids optimal et on a donc une perte de précision de l'estimateur MGPP double par rapport à l'estimateur MGPP simple optimal.

L'intérêt de la MGPP double, comme utilisée par La Poste, est de diminuer le nombre de liens à observer pour calculer les poids standardisés. Cependant, cet avantage se limite aux liens de type TpU-TpU. Pour ce type de liens, on s'intéresse à deux méthodes de standardisation :

— La standardisation double, $\sum_{i \in U_F} \hat{\theta}_{ij}^{FM} = 1$ et $\sum_{j \in U_M} \hat{\theta}_{jk}^{MT} = 1$.

— D'autres types de standardisation globale, $\sum_{i \in U_F} \sum_{j \in U_M} \hat{\theta}_{ij}^{FM} \hat{\theta}_{jk}^{MT} = 1$, moins contrai-
gnants.

Dans l'article Medous et al. (2023a), plusieurs méthodes de standardisation sont comparées : la standardisation de la MGPP simple (voir section 1.1.4), la standardisation double, et un cas particulier de standardisation globale, où $\theta_{ik} = \sum_{j \in U_M} \theta_{ij}^{FM} \theta_{jk}^{MT}$, $i \in U_F$, $k \in U$ est divisé par le total $\sum_{i \in U_F} \theta_{ik}$. Soit $N_{Fk}$ le nombre de liens entre $U_F$ et $k \in U$, $N_{Mk}$ le nombre de liens entre $U_M$ et $k \in U$ et $N_{Fj}$ le nombre de liens entre $U_F$ et $j \in U_M$. Dans le cas de La Poste, $N_{Fk}$ est le nombre d'adresses dans la tournée $k$, $N_{Mk}$ le nombre de cases dans la tournée $k$ et $N_{Fj}$ le nombre d'adresses dans la case $j$. Si les liens sont de type TpU-TpU, $N_{Fk} = \sum_{j \text{ lié à } k} N_{Fj}$, c'est-à-dire que le nombre d'adresses dans une tournée est égal à la somme des adresses contenues dans toutes les cases de la tournée. Il faut observer $N_{Fk}$ liens pour la standardisation de la MGPP simple, $N_{Fk} + N_{Mk}$ liens pour la standardisation globale considérée et $N_{Fj} + N_{Mk}$ pour la standardisation double (voir l'exemple de la figure 1.26). En pratique, si $N_{Fj} + N_{Mk} < N_{Fk}$, alors la standardisation double permet un gain en terme de liens à observer. En particulier, si $N_{Fk} = N_{Fj}N_{Mk}$ et $N_{Fj} > 2$, $N_{Mk} > 2$, alors $N_{Fj} + N_{Mk} < N_{Fk}$ et la MGPP double permet de réduire le nombre de liens à observer par rapport à la MGPP simple. Dans le cas de La Poste, la condition $N_{Fk} = N_{Fj}N_{Mk}$ est satisfaite si toutes les cases contiennent le même nombre d'adresses, ce qui n'est pas
Type de standardisation

MGPP simple :

$$
\sum_{i \in U_F} \tilde{\theta}_{ik} = 1
$$

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

nombre de liens observés: 10

globale :

$$
\sum_{i \in U_F} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} = 1
$$

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

nombre de liens observés: 12

$$
\begin{array}{ll}
\text{double} : & \sum_{i \in U_F} \tilde{\theta}_{ij}^{FM} = 1 \\
\text{et} & \sum_{j \in U_M} \tilde{\theta}_{jk}^{MT} = 1
\end{array}
$$

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

nombre de liens observés: 7

<caption>FIGURE 1.26 – Comparaison des méthodes de standardisation.</caption>

réaliste. Par contre, le nombre d'adresses par case plus le nombre de cases par tournée vaut en moyenne 60 alors que le nombre d'adresses par tournée vaut en moyenne 500. La condition $N_{Fj} + N_{Mk} < N_{Fk}$ est donc satisfaite pour un large nombre de tournées et de cases, ce qui permet un gain important en terme de nombre de liens à observer.

On a vu précédemment que la différence entre la variance de l'estimateur **MGPP** simple quelconque et la variance de l'estimateur **MGPP** simple optimal du total dépend, via la standardisation, de la répartition des liens entre les populations. Pour mieux comprendre l'impact de la répartition des liens sur la précision de la **MGPP** utilisée à La Poste, deux études par simulation ont été faites. Dans chaque simulation, la population cible est la population des tournées, la population intermédiaire celle des cases et la population Frame celle des adresses.

La première étude montre que l'augmentation de variance de la **MGPP** double compa-
rée à la **MGPP** simple optimale dans le cas de La Poste est négligeable quand le nombre
d'adresses par case est constant, indépendamment du nombre de cases par tournée. Si le
nombre d'adresses par case varie de façon importante selon les cases, la perte de précision
de la **MGPP** double dépend du nombre de cases par tournée. La perte de précision est alors
plus importante si le nombre de cases varie selon les tournées.

La deuxième étude cherche à reproduire la situation de La Poste pour évaluer la perte
de précision de la MGPP double par rapport à la MGPP simple optimale, mais aussi de la
MGPP simple optimale par rapport à l'estimateur direct. Les populations sont générées de
manière à reproduire la situation observée sur les historiques de La Poste. Cette simulation montre que la MGPP optimale est un peu moins précise que l'estimateur direct et que la MGPP double est bien moins précise que la MGPP optimale. L'utilisation d'une MGPP double non optimale (voir section 1.2.1) par La Poste est donc responsable de la détérioration de la variance des estimateurs du trafic postal observée depuis 2015 (voir section 1.2.1). L'utilisation de la MGPP simple optimale dans le cas de La Poste pourrait permettre de conserver une variance proche de celle de l'estimateur direct. Ce dernier point fait l'objet de l'article suivant.

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

Dans l'article Medous et al. (2023a), précédemment résumé, j'introduis le concept de MGPP double et je discute des gains en terme de liens observés quand les liens sont de type TpU. Je note aussi que la MGPP double utilisée par La Poste induit une perte de précision par rapport à la MGPP optimale.

Dans le chapitre 3, je me pose la question de l'existence de poids optimaux pour la MGPP double avec double standardisation, c'est à dire des poids standardisés appelés $\tilde{\theta}_{ij}^{FM, opt}$ et $\tilde{\theta}_{jk}^{MT, opt}$, $i \in U_F$, $j \in U_M$, $k \in U$, tels que la somme $\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{jk}^{MT, opt}$ satisfait l'équation (1.10) :

$$ \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \sum_{k \in U} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM, opt} \tilde{\theta}_{jk}^{MT, opt} y_k \sum_{k' \in U} \sum_{j' \in U_M} \tilde{\theta}_{i'j'}^{FM, opt} \tilde{\theta}_{j'k'}^{MT, opt} y_{k'} \\ \le \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \sum_{k \in U} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} y_k \sum_{k' \in U} \sum_{j' \in U_M} \tilde{\theta}_{i'j'}^{FM} \tilde{\theta}_{j'k'}^{MT} y_{k'} $$

pour toute variable d'intérêt $y_k$, $k \in U$, et tout poids standardisé $\tilde{\theta}_{ij}^{FM}$ et $\tilde{\theta}_{jk}^{MT}$, $i \in U_F$, $j \in U_M$, $k \in U$. Ce problème est plus complexe que celui de l'existence de poids optimaux pour la MGPP simple, car il faut étudier deux ensembles de liens. Par souci de simplicité, on désignera la MGPP double donnée par l'équation (1.14), avec double standardisation, par MGPP double dans la suite de ce chapitre et dans le chapitre 3.

Je montre que des poids optimaux pour la MGPP double existent pour des liens de type TpU-TpU ("Tous-pour-Un-Tous-pour-Un"), pour des plans de sondage satisfaisant les conditions détaillées ci-dessous. Ces conditions sont plus fortes que la $\Delta$-propriété pour U telle que définie dans Medous et al. (2023a). La première condition est la $\Delta$-propriété entre $U_F$ et $U_M$ et la seconde, appelée $c_k$-condition, porte sur les liens entre $U_M$ et $U$.

La $\Delta$-propriété définie entre $U_F$ et $U$ dans Medous et al. (2023a) se définit similairement entre les populations $U_F$ et $U_M$. Soit $N_F$ la taille de la population Frame $U_F$ et, pour $j \in U_M$, $U_{Fj}$ la sous-population de $U_F$ de taille $N_{Fj}$ composée des individus reliés à j. On suppose
dans ce qui suit la population $U_F$ ordonnée selon les sous populations $U_{Fj}$, $j \in U_M$, et on note
$\Delta_{jj'}^M = (\Delta_{ii'})_{i \in U_{Fj}, i' \in U_{Fj'}}$ la sous-matrice de $\Delta$ de taille $N_{Fj} \times N_{Fj'}$ composée des éléments
en position $i, i'$ si $i$ (resp. $i'$) est lié à $j$ (resp $j'$). Soit $1_j^M$ le vecteur de taille $N_{Fj}$ dont les
composantes sont égales à 1. Si les liens entre $U_F$ et $U_M$ sont de type TpU, on définit la
$\Delta$-propriété pour $U_M$ par :

- pour tout $j \in U_M$, $\Delta_j^M$ est inversible,

- pour $j \neq j' \in U_M$,

$$ \Delta_{j,j' \neq j}^{M} = c_{jj'} \mathbf{1}_{j}^{M} \mathbf{1}_{j'}^{Mt} \quad \text{avec } c_{jj'} \text{ ne dépendant ni de } i \text{ ni de } i'. \quad (1.15) $$

La $\Delta$-propriété pour $U_M$ permet d'obtenir la décomposition de $\Delta$ présentée dans la figure 1.27, avec $\Delta_{kk'}^T$, $k \in U$ les sous-matrices de $\Delta$ utilisées pour la $\Delta$-propriété entre $U_F$ et $U$ (voir chapitre 2). Cette propriété n'implique pas la $\Delta$-propriété pour $U$ et on rajoute une condition sur les liens entre $U_M$ et $U$ pour obtenir des poids optimaux.

Pour $k \in U$, on définit par $U_{Mk}$ la sous-population de $U_M$ de taille $N_{Mk}$ contenant tous les individus $j \in U_M$ liés à $k$. On suppose les individus de $U_M$ ordonnés selon ces populations (voir la figure 1.27). Soient les sous-populations $U_{Fk}$, $k \in U$, de $U_F$, définies dans Medous et al. (2023a), avec $U_{Fk}$ contenant les individus de $U_F$ reliés à $k \in U$. Si les liens sont de type TpU-TpU, $U_F$ est ordonnée selon les sous-population $U_{Fj}$, $j \in U_M$, et $U_M$ est ordonnée selon les sous-populations $U_{Mk}$, $k \in U$, alors $U_F$ est ordonnée selon les sous-populations $U_{Fk}$, $k \in U$.

Si les liens sont de type TpU-TpU, on dit qu'un plan de sondage satisfait la $c_k$-condition s'il satisfait la $\Delta$-propriété pour $U_M$ et que, pour tout $j$ et $j' \in U_M$ liés respectivement à $k$ et $k' \in U$ tels que $j \neq j'$, la constante $c_{jj'}$ donnée par la $\Delta$-propriété pour $U_M$ vérifie :

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

avec $c_{kk'}$ ne dépendant ni de $j$ ni de $j'$ et $c_{kk} = c_k$. La $c_k$-condition est vérifiée pour les plans de Poisson et SASSR, ainsi que, sous certaines conditions (détaillées dans le chapitre 3), pour le STASSR.

Si la $c_k$-condition est vérifiée, alors je montre qu'il existe des poids minimisant la variance de l'estimateur MGPP double pour toute variable d'intérêt et que ces poids optimaux prennent la valeur suivante :
<figure><img src="image_34.png" /><figcaption>FIGURE 1.27 – Décomposition de la matrice Δ dans le cas TpU-TpU.</figcaption></figure>

$$ (\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} $$

$$ \tilde{\theta}_{jk}^{MT, opt} = \begin{cases} \frac{\mathbf{1}_j^{Mt}(\Delta_j^M)^{-1} \mathbf{1}_j^M}{\sum_{j' \in U_{Mk}} \mathbf{1}_{j'}^{Mt}(\Delta_{j'}^M)^{-1} \mathbf{1}_{j'}^M} / (1 - c_k \mathbf{1}_j^{Mt}(\Delta_j^M)^{-1} \mathbf{1}_j^M) & \text{si } j \text{ est lié à } k, \\ 0 & \text{sinon} \end{cases} \quad (1.16) $$

où la valeur de $c_k$ est donnée par la $c_k$-condition et $\tilde{\theta}_{ij}^{FM, opt}$, $i \in U_F$, $j \in U_M$, sont les poids optimaux pour la MGPP simple entre $U_F$ et $U_M$, donnés par l'équation (1.12).

La différence entre la variance de l'estimateur MGPP double utilisant des poids $\tilde{\theta}_{ij}^{FM}$ et
$\tilde{\theta}_{jk}^{MT}$, $i \in U_F$, $j \in U_M$, $k \in U$, quelconques et la variance de l'estimateur MGPP double utilisant les poids optimaux définis par l'équation (1.16) est obtenue en remplaçant les poids $\tilde{\theta}_{ik}^{opt}$ par $\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}$ dans l'équation (1.13).

Cependant, si on compare un estimateur MGPP $\hat{t}_{y2}^{FM, opt}$ obtenu en utilisant les poids optimaux $\tilde{\theta}_{ij}^{FM, opt}$ donnés dans (1.16) et des poids quelconques $\tilde{\theta}_{jk}^{MT}$ avec l'estimateur MGPP double optimal, alors la différence de variances se simplifie comme suit :

$$ \mathrm{Var}(\hat{t}_{y2}^{FM, opt}) - \mathrm{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})^2 \left( \frac{1}{1_j^{Mt}(\Delta_j^M)^{-1} 1_j^M} - c_k \right). \quad (1.17) $$

Les poids optimaux pour la MGPP double donnés dans (1.16) permettent d'obtenir la même précision que l'estimateur MGPP simple optimal tout en réduisant le nombre de liens observés. Cependant, ces poids ne sont calculables que si l'information $1_j^{Mt}(\Delta_j^M)^{-1} 1_j^M / (1 - c_k 1_j^{Mt}(\Delta_j^M)^{-1} 1_j^M)$ est disponible pour $j \in U_M$. Pour le plan SASSR, les poids optimaux prennent une forme très simple :

$$ (\tilde{\theta}_{ij}^{FM, opt})_{i \in U_{Fj}} = 1/N_{Fj} \text{ si } i \text{ et } j \text{ sont liés, 0 sinon, et} $$

$$ \tilde{\theta}_{jk}^{MT, opt} = N_{Fj}/N_{Fk} \text{ si } j \text{ et } k \text{ sont liés, 0 sinon,} $$

avec $N_{Fj}$ le nombre d'individus $i \in U_F$ liés à $j \in U_M$ et $N_{Fk}$ le nombre d'individus $j \in U_M$ liés à $k \in U$. Si les liens sont de type TpU-TpU, $N_{Fk} = \sum_{j \in U_{Mk}} N_{Fj}$. Il faut donc connaître $N_{Fj}$ pour tout $j \in U_M$ lié à $k$ pour calculer $\tilde{\theta}_{jk}^{MT, opt}$.

Pour La Poste, $N_{Fj}$ est le nombre d'adresses dans la case $j$ et $N_{Fk}$ le nombre d'adresses dans la tournée $k$. Pour calculer les poids $\tilde{\theta}_{jk}^{MT, opt}$, il faut connaître le nombre $N_{Fj}$ d'adresses contenues dans la case $j$ pour toutes les cases $j \in U_M$ liées à la tournée $k$. Cette information n'est pas disponible et il faut, pour l'obtenir, compter les adresses contenues dans les cases liées à la tournée $k$. Cela revient à compter toutes les adresses desservies par la tournée $k$ et La Poste fait donc face au même problème qu'avec la MGPP simple.

Pour pallier ce problème, je propose des poids alternatifs pour la MGPP double, qui permettent de limiter la perte de précision des estimateurs comparés aux estimateurs MGPP double optimaux. En pratique, les poids optimaux $\tilde{\theta}_{ij}^{FM, opt}$, $i \in U_F$, $j \in U_M$ sont déjà calculés par La Poste dans SYCI 2, je me concentre donc sur des poids alternatifs aux poids optimaux $\tilde{\theta}_{jk}^{MT, opt}$.

Soit une variable auxiliaire $x$ définie sur $U_M$, dont la valeur $x_j$ est disponible pour tout
$j \in U_M$. On définit les poids alternatifs suivants pour tout $j \in U_M$ et $k \in U$ :

$$ \tilde{\theta}_{jk}^{MT,alt} = \begin{cases} \frac{x_j}{\sum_{j' \in U_{Mk}} x_{j'}} & \text{si } j \text{ lié à } k, \\ 0 & \text{sinon.} \end{cases} $$

En me basant sur la différence de variances donnée par (1.17), je propose deux façons de sélectionner la variable auxiliaire $x$ utilisée pour calculer les poids alternatifs $\tilde{\theta}_{jk}^{MT,alt}$.

La première méthode consiste à choisir $x$ corrélée au poids optimal non standardisé $(\theta_{jk}^{MT, opt})_{j \in U_M, k \in U}$. Dans le cas de La Poste, cela revient à sélectionner la variable $x$ la plus corrélée au nombre d'adresses par case. La Poste observe le nombre d'adresses par case et d'autres informations sur les cases (voir section 1.2.2) pour les cases indirectement échantillonnées dans SYCI 2. La corrélation est donc facilement calculable à partir de l'échantillon ou, dans le cas de La Poste, à partir d'un historique de données. Cependant, cette méthode ne prend pas en compte l'impact de $y$ sur la différence de variances et peut donc donner de mauvais résultats si $y$ prend de grandes valeurs (voir section 3.5).

Les liens étant de type TpU, il n'y a qu'une unité $k \in U$ liée à $j \in U_M$ et il est possible de calculer les valeurs $\tilde{y}_j = y_k \tilde{\theta}_{jk}^{MT, opt}$ pour les unités de $U_M$. La seconde méthode proposée consiste à choisir $x$ corrélée à la variable $\tilde{y}$ prenant les valeurs $\tilde{y}_j$ pour $j \in U_M$. Cette méthode permet de prendre en compte l'impact de $y$, mais la valeur $\tilde{\theta}_{jk}^{MT, opt}$ étant inconnue, l'hypothèse que la variable $x$ est corrélée à $\tilde{y}$ est invérifiable en pratique. De plus, la variable $x$ choisie dépend de la variable d'intérêt et il faut recalculer les poids alternatifs si la variable sélectionnée change selon les variables d'intérêt.

J'appuie cette discussion sur des simulations basées sur les données de La Poste (voir section 3.5.3). Cette étude par Monte-Carlo montre que la première méthode donne des estimateurs plus précis que la deuxième méthode quand les corrélations sont élevées, alors que la deuxième méthode donne des estimateurs plus précis quand les corrélations sont faibles. Elles montrent aussi que, dans le cas de La Poste, le seuil de corrélation pour que la deuxième méthode donne des estimateurs plus précis que la MGPP double actuelle est faible, alors qu'il est plus élevé avec la première méthode. Cette étude permet une meilleure compréhension du problème de réduction de la variance pour la MGPP double, mais il reste difficile d'améliorer la précision des estimateurs à La Poste à l'heure actuelle. Certaines bases de données (dont RAO, voir section 1.2.2) contiennent des variables auxiliaires qui semblent fortement corrélées au nombre d'adresses par case, comme un nombre d'adresses par case théorique, mais qui ne sont pas encore exploitables.
### 1.3.2 Intégration statistique de données

Une part importante de mes travaux de thèse a porté sur l'utilisation des bases de données massives dont dispose La Poste, avec en priorité la base TAE (voir section 1.2.2). Dans la section précédente, on a montré que des données auxiliaires peuvent être utilisées dans le calculs des poids de la MGPP double. Les données auxiliaires sont alors utilisées pour améliorer la méthodologie actuelle sans la modifier. Cependant la réduction de variance apportée par cette méthode est limitée, car on peut, au mieux, obtenir l'estimateur MGPP double optimal.

Dans cette section, je propose une nouvelle méthodologie basée sur l'utilisation des données disponibles dans TAE. Les résultats présentés seront publiés dans l'article "QR Prediction for Statistical Data Integration" (Survey Methodology, à paraître).

#### Chapitre 4 : QR Prediction for Statistical Data Integration

La littérature sur l'utilisation de bases de données massives en sondage considère le cas où les valeurs des variables d'intérêt sont connues pour tous les individus des bases de données non probabilistes (Kim and Tam, 2021; Beaumont, 2020). Or, les données de La Poste, TAE (voir Section 1.2.2), sont collectées par les machine de tri. Certaines variables d'intérêt considérées par La Poste, comme le type ou le montant d'affranchissement, ne sont pas reconnues par les machine et leurs valeurs pour les individus de TAE ne sont donc pas connues. Dans cet article, je commence par étudier certains estimateurs proposés par Kim and Tam (2021), puis je propose un estimateur du total $T_y$ adapté aux données postales inspirés de l'estimateur QR introduit par Wright (1983) et Särndal and Wright (1984).

Soit $s_P$ un échantillon probabiliste tiré avec un plan de sondage $p(.)$ et $s_{NP}$ un échantillon non-probabiliste de taille $N_{NP}$. On note $\pi_k, k \in U$ (resp. $\pi_{kk'}, k, k' \in U$) les probabilités d'inclusion dans $s_P$ de premier (resp. second) ordre et $\delta_k$ l'indicatrice d'inclusion dans $s_{NP}$, où $\delta_k = 1$ si $k$ appartient à $s_{NP}$, 0 sinon. On suppose que la valeur de $\delta_k$ est connue pour tout individu $k$ appartenant à $s_P$.

On suppose dans un premier temps que la variable d'intérêt $y$ est connue pour tout individu de $s_P$ et de $s_{NP}$. On note $T_{NP} = \sum_{k \in s_{NP}} y_k = \sum_{k \in U} \delta_k y_k$ le total de $y$ sur $s_{NP}$ et $T_C = T_y - T_{NP}$ le total de $y$ sur $U - s_{NP}$. Si la taille de la population $N$ est connue, Kim and Tam (2021) proposent l'estimateur suivant :

$$ \hat{T}_{PDI} = T_{NP} + \hat{T}_{C}^{(Ha)}, \qquad (1.18) $$
avec $\hat{T}_C^{(\text{Ha})}$ l'estimateur de Hájek du total $T_C$ (voir section 1.1.3) :

$$ \hat{T}_C^{(\text{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)}. $$

Je montre que la différence entre la variance de l'estimateur de **HT** du total $\hat{T}_{HT,y}$ obtenu à partir de l'échantillon $s_P$ (voir section 1.1.3) et la variance asymptotique de $\hat{T}_{PDI}$ est donnée par :

$$ \text{Var}(\hat{T}_{HT,y}) - \text{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), $$

où $\bar{Y}_U = \frac{1}{N} \sum_{k \in U} y_k$ est la moyenne de $y$ sur $U$, et $\bar{Y}_C = \frac{1}{N - N_{NP}} \sum_{k \in U} (1 - \delta_k)y_k$ est la moyenne de $y$ sur $U - s_{NP}$.

L'estimateur $\hat{T}_{PDI}$ est asymptotiquement plus précis, en terme de variance, que l'estimateur de Horvitz-Thompson $\hat{T}_{HT,y}$. La Poste souhaite utiliser l'estimateur $\hat{T}_{PDI}$, en utilisant les échantillons probabilistes tirés dans le cadre de **SYCI 2** et **TAE** comme échantillon non-probabiliste (voir 1.2.2), afin de réduire la variance de ses estimateurs actuels et limiter l'impact de la réduction future des tailles d'échantillons (voir section 1.2.1). Cependant, La Poste ne dispose pas la valeur de $y$ pour les individus de **Traitement Automatisé de l'Enveloppe (TAE)**. La Poste ne peut donc pas utiliser $\hat{T}_{PDI}$ tel quel.

Dans la deuxième partie de Medous et al. (2023b), on suppose que les valeurs de $y$ sont connues pour les individus de $s_P$ mais inconnues pour les individus de $s_{NP}$. On suppose aussi que l'intersection entre $s_P$ et $s_{NP}$ est non vide, ce qui est le cas à La Poste car la base de données **TAE** contient 90% des lettres distribuées (voir section 1.2.2). Dans cette partie, je cherche à adapter l'estimateur $\hat{T}_{PDI}$ aux données de La Poste.

On suppose qu'un vecteur $\boldsymbol{x}_k = (X_{k1}, \dots, X_{kp})^\top$ de variables auxiliaires est disponible pour tout individu $k \in s_{NP}$ et que les valeurs de $\delta_k$ et $\delta_k \boldsymbol{x}_k$ sont connues pour tout individu $k \in s_P$.

L'idée présentée dans Medous et al. (2023b) est de remplacer les valeurs $y_k$, $k \in s_{NP}$ dans l'équation (1.18) par des prédictions $\hat{y}_k$, obtenues en modélisant $y$ en fonction des variables auxiliaires sur l'intersection $s_P \cap s_{NP}$. On obtient alors un prédicteur $\hat{T}_{NP}$ du total $T_{NP}$.

On suppose le modèle suivant :

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

où les erreurs $\varepsilon_k$ sont indépendantes, d'espérance $E_m(\varepsilon_k) = 0$ et de variance $\text{Var}_m(\varepsilon_k)$ propor-
tionnelle à $\nu(\boldsymbol{x}_k) = v_k$ pour une constante positive $v_k$. L'indice $m$ indique que l'on considère l'espérance et la variance sous le modèle 1.19.

Soit $q_k, k \in U$ des constantes positives connues. On définit le prédicteur de $y_k, k \in U$ suivant :

$$ \hat{y}_k = \boldsymbol{x}_k^\top \hat{\boldsymbol{\beta}} $$

avec

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

où la matrice $\sum_{k \in s_P} q_k \delta_k \boldsymbol{x}_k \boldsymbol{x}_k^\top$ est supposée inversible pour tout échantillon $s_P$. Je propose dans Medous et al. (2023b) d'utiliser un prédicteur de type QR, introduit par 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 \boldsymbol{x}_k^\top \hat{\boldsymbol{\beta}} + \sum_{k \in s_P} r_k \delta_k (y_k - \boldsymbol{x}_k^\top \hat{\boldsymbol{\beta}}), \end{aligned} \quad (1.21) $$

où $r_k, k \in U$ sont des constantes positives connues.

Ce prédicteur est sans biais sous le modèle 1.19 (voir chapitre 4). En pratique, on préfère détailler le biais du prédicteur QR sous le plan.

Afin de pouvoir étudier le comportement du prédicteur $\hat{T}_{NP}^{(QR)}$ sous l'inférence du plan de sondage, on suppose qu'il existe un vecteur $\boldsymbol{\lambda} \in \mathbf{R}^p$ tel que

$$ (QR): \quad 1 - \pi_k r_k = \pi_k q_k \boldsymbol{x}_k^\top \boldsymbol{\lambda} \quad \text{pour tout} \quad k \in s_{NP}. \quad (1.22) $$

Cette condition, appelée (QR) condition dans Medous et al. (2023b), permet de simplifier l'expression du prédicteur QR comme suit :

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

avec

$$ \hat{T}_{NP}^{(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}}) \quad (1.23) $$

où $d_k = 1/\pi_k$.

Wright (1983) et Särndal and Wright (1984) ont montré que, sous certaines conditions, le prédicteur QR défini par (1.23) est asymptotiquement sans biais. Les conditions nécessaires,
ainsi que les preuves, sont détaillées dans le chapitre 4.

Le total $T_y$ est estimé par :

$$ \hat{T}_y^{(QR)} = \hat{T}_{NP}^{(QR)} + \hat{T}_C^{(Ha)}. \qquad (1.24) $$

Dans le chapitre 4, je montre que $\hat{T}_y^{(QR)}$ a des propriétés de biais similaires à celles de $\hat{T}_{NP}^{(QR)}$.

Si la condition (QR) est vérifiée, et sous des hypothèses sur le plan de sondage et les variables y et $\boldsymbol{x}$, l'erreur d'estimation peut être approchée par :

$$ \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) = o_p\left(\frac{1}{\sqrt{n}}\right), $$

avec

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

et

$$ \hat{\beta} = \left( \sum_{k \in U} \pi_k q_k \delta_k \boldsymbol{x}_k \boldsymbol{x}_k^\top \right)^{-1} \sum_{k \in U} \pi_k q_k \delta_k \boldsymbol{x}_k y_k. $$

Cette approximation permet d'obtenir la variance approchée de l'estimateur QR :

$$ \mathrm{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). $$

Si $\pi_{kk'} > 0$ pour tout $k, k' \in U$, un estimateur de cette variance est donné par :

$$ \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), $$

avec $\hat{E}_k = \delta_k(y_k - \boldsymbol{x}_k^\top \hat{\beta})$ et $\hat{e}_k = (1 - \delta_k)(y_k - \sum_{k' \in s_p} d_{k'}(1 - \delta_{k'}))y_{k'}/(N - N_{NP})$.

Dans le chapitre 4, on discute du choix des constantes $q_k$ et $r_k$ et on propose trois options qui vérifient la condition (QR) donnée par (1.22) :

$$ q_k = d_k v_k^{-1} \text{ et } r_k = d_k, $$

$$ q_k = v_k^{-1} \text{ et } r_k = 1, $$

$$ q_k = (d_k - 1) v_k^{-1} \text{ et } r_k = 1. $$

Les deux premières options permettent d'obtenir respectivement l'estimateur Model-Assisted et Model-Based (Särndal et al., 1989). La troisième option permet d'obtenir l'estimateur
cosmétique du total $T_{NP}$, introduit par Särndal and Wright (1984) :

$$ \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})}), $$

avec $\hat{y}_k^{(\text{Cos})} = x_k^\top \hat{\beta}^{(\text{Cos})}$ et

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

Une étude sur populations simulées permet de mettre en avant les avantages et inconvénients de chaque option et de comparer leur efficacité (en terme de biais et variance) par rapport à l'estimateur de **HT** et l'estimateur $\hat{T}_{PDI}$ du total $T_y$.

Une deuxième étude utilisant les données historiques de La Poste montre des résultats similaires. Dans le cas de La Poste, l'utilisation d'un estimateur cosmétique du total $T_{NP}$ semble la plus pertinente. Cependant, on note que l'utilisation de l'estimateur de Hájek du total $T_C$ n'est pas la plus adaptée aux données de La Poste, d'une part à cause de sa sensibilité à la sélection de $s_P$ (voir chapitre 4) et d'autre part parce que la taille $N$ de la population à La Poste est inconnue. L'utilisation d'un estimateur de **HT** du total $T_C$ est plus judicieuse pour La Poste et est actuellement à l'étude.
Sondage Indirect et Méthode
Généralisée de Partage des Poids
<figure><img src="image_35.png" /></figure>
# Chapitre 2

## Many-to-One indirect sampling with
application to the French postal traffic
estimation<fnref n="1" />, modifié pour la thèse

### Abstract

In social and economic surveys, it can be difficult to directly reach units of the target population, and indirect sampling is often advocated to solve this issue. In indirect sampling, the sample is drawn from a frame population that is linked to the target population, and estimation of target population parameters is typically achieved through the **Generalized Weight Share Method (GWSM)**. This method provides a weight, for every unit of the target population, that depends on the one hand, on the sampling weights in the frame population and, on the other hand, on the link weights between the frame population and the target population. In the present study, we focus on the situation in which the units from the frame population are linked to one and only one unit from the target population (Many-to-One case). This situation is encountered at the French postal service where addresses are sampled instead of postman rounds. We aim at understanding of the impact of the link weights on the efficiency of the **GWSM** estimators. We derive variance expressions and optimality results for a large class of sampling designs. Moreover, we note that the Many-to-One case can lead to too many links to observe. We alleviate the problem by introducing an intermediate population and double indirect sampling. The question of the loss of precision in this situation is discussed in detail through theoretical results and simulations. These findings help to explain the loss of precision of double **GWSM** estimators observed recently at the French postal service.

<footnote n="1">1. Article de E.Medous, C.Goga, A. Ruiz-Gazen, J.F Beaumont, A. Dessertaine et P. Puech, 2023, publié
dans *The Annals of Applied Statistics*, 17(1):838-859</footnote>
**Keywords** : Generalized Weight Share Method, Optimal link weights, Stratified sam-
pling, Variance estimation.

**2.1 Introduction**

In France, at the postal service (La Poste), only part of the postal traffic goes through an automatized processing. The monthly postal traffic is unknown and is estimated through probability-based surveys. For many years, La Poste has drawn samples directly in the population of postman rounds, which is considered to be the population of interest or the target population. Since 2008, the organization of postman rounds has changed and is no longer stable over time. Sampling directly in the target population has become impossible, and the sampling design has evolved to an indirect sampling design where the frame population is the population of postal addresses.

Indirect sampling has been extensively studied in the survey sampling literature; see
Deville and Lavallée (2006), Lavallée (2007), Kiesl (2016) and Haziza and Beaumont (2017)
for a general theory and reviews with many references inside. As described in Kiesl (2016),
indirect sampling has many applications, including household panels : Kalton and Brick
(1995); Rendtel and Harms (2009) and hard-to-reach populations : Deville and Maumy-
Bertrand (2006) for tourism and De Vitiis et al. (2014) for the homeless population. The use
of indirect sampling at La Poste, has been introduced in Dessertaine and Fluteaux (2004) and
Lardin-Puech (2014). A useful estimation method in the context of indirect sampling is the
**Generalized Weight Share Method (GWSM)**, as detailed in Deville and Lavallée (2006). It
consists of using the links that relate the frame and the target populations, and considering a
total over the target population to be a total over the frame population. The use of standard
methods, such as the Horvitz-Thompson estimator, is then possible and leads to the **GWSM**
estimator. Kalton and Brick (1995), Deville and Lavallée (2006) and Kiesl (2016), among
others, studied in detail the properties of the **GWSM** estimator, and in particular the question
of the impact of the link structure and link weights on its variance. Optimality, in the sense
of variance minimization with respect to the link weights for unbiased **GWSM** estimators,
is discussed at length in Deville and Lavallée (2006). The conclusion is that optimal **GWSM**
estimators, that do not depend on the variable of interest, cannot be derived for general link
structure. In the present study, we propose to focus on a particular link structure described
below that is of interest at La Poste, and go further into the understanding of indirect
sampling.

At the French postal service, every postal address is linked to only one postman round
for a given day. This link structure is of a particular type, called Many-to-One (MtO), where each unit in the frame population is linked to one and only one unit in the target population. This situation is also encountered in households surveys where individuals are sampled instead of households. This link structure is studied in detail in the present paper. We derive the optimal GWSM estimator that minimizes the variance among the unbiased GWSM estimators, for a large class of indirect MtO sampling designs. This class includes Poisson sampling, simple random sampling without replacement, and stratified designs, including the design implemented at La Poste. Moreover, we derive a simple formula to evaluate the increase in variance when using a non-optimal GWSM estimator compared to the optimal one.

The weight share method is simple but requires that the links between the indirectly sampled units in the target population, and the frame population, are known. The problem faced by the French postal service with MtO links, is that every unit in the target population is linked to a very large number of units in the sampling frame. At La Poste, all addresses delivered during a sampled postman round must be known. On average, there are approximately 500 addresses per postman round, and it is not possible to enumerate all of the addresses in the morning, before the departure of the postman.

To get around this problem, La Poste has set up a double indirect sampling design, using the outgoing mail sorting boxes as an intermediate population. This method is much faster than simple indirect sampling. Only the addresses of the sampled boxes and the boxes of the sampled rounds are to be observed, which is approximately 60 items on average, compared to the 500 for simple indirect sampling. Given the situation at La Poste, double indirect sampling is an alternative to a time consuming simple indirect sampling design. This alternative is necessary to be able to collect the data. However, using this method, La Poste observed a deterioration in the precision of the estimators. The goal of the present paper is to understand the loss of precision observed at La Poste, but also to give guidelines for the implementation of an efficient double indirect sampling design.

In Section 2.2 of the present paper, we consider a large class of indirect MtO sampling designs. We derive the optimal GWSM estimator and give a simple expression of the difference in the variances between the optimal GWSM estimator and any non-optimal unbiased GWSM estimator. For Poisson sampling, we also prove that the optimal GWSM estimator is less precise than the direct estimator. The result on optimal GWSM estimator is used in Section 2.3, where we introduce and compare double indirect MtO sampling with simple indirect MtO sampling with optimal link weights. In the same section, we detail situations where there is a gain, in terms of the smaller number of links to observe, from using double indirect sampling compared to simple indirect sampling. In Section 2.4, we define several
setups and illustrate, through a Monte Carlo study, the impact of double indirect sampling
on the precision of the estimators. Depending on the link structure, we observe that there
could be no loss at all, or, on the contrary, an enormous loss of precision. In Section 2.5,
we also give numerical results in a context similar to La Poste. These results allow us to
explain the loss of precision observed at La Poste, when using a double indirect sampling
design compared to a simple direct sampling design. Section 2.6 concludes the paper while
the proofs are given in the Appendix.

## 2.2 Indirect sampling

2.2.1 GWSM

In some surveys, it is not possible to sample directly from the target population $U_T$.
However, a sampling frame can exist for a population $U_F$, that is related to $U_T$ in such a
way that any unit in $U_T$ is linked to at least one unit in $U_F$. Indirect sampling refers to
selecting a sample $s_F$ from $U_F$ by using standard selection methods and derive estimators
for parameters defined on $U_T$. In the case of the La Poste survey, the population of interest
is made of postman rounds on a given day in France, but no sampling frame for the rounds
exists. A sampling frame for postal addresses is however available, and each postman round
contains at least one address.

Let us denote by $N_T$ (resp. $N_F$) the size of $U_T$ (resp. $U_F$) and by $l_{ik}$ the link between
$i \in U_F$ and $k \in U_T$, with $l_{ik} = 1$ if the units $i$ and $k$ are linked, and $l_{ik} = 0$ otherwise. Units
from $U_F$ can be linked in several ways to units from $U_T$ (Deville and Lavallée, 2006). We
can have “Many-to-One” (MtO) links as on the left panel of Figure 2.1, namely each unit
from the frame population $U_F$ is linked to only one unit from the target population $U_T$. We
can have “One-to-Many” links as on the middle panel of Figure 2.1, namely each unit from
$U_T$ is linked to only one unit from $U_F$. Finally, we can have “Many-to-Many” (MtM) links
as on the right panel of Figure 2.1, with units from $U_F$ linked to several units in $U_T$ and
reciprocally. In the La Poste survey, an address almost always belongs to only one round and
the links are MtO. Following Deville and Lavallée (2006), we start by making the assumption
that the links between $U_F$ and $U_T$ can be observed for every unit in $U_F$ and every unit in
$U_T$.

Let $y$ be the variable of interest measured on $U_T$, and let $y_k$ be its value for the $k$-th
unit in $U_T$. We are interested in estimating $t_y = \sum_{k \in U_T} y_k$, the total of $y$ over $U_T$. A sample
$s_F$ is drawn from $U_F$ according to a sampling design $p_F(\cdot)$. We can associate to $s_F$ the
vector $(I_1, \dots, I_{N_F})'$ where $I_i$ is the sample membership indicator of the individual $i$ from
<figure><img src="image_36.png" /><figcaption>FIGURE 2.1 – The different types of links between $U_F$ and $U_T$.</figcaption></figure>

$U_F$ defined as $I_i = 1$ if $i$ is selected and $I_i = 0$ otherwise. We denote by $\pi_i = \Pr(i \in s_F)$ the first-order inclusion probability of unit $i$ and by $\pi_{ii'} = \Pr(i, i' \in s_F)$ the second-order inclusion probability of units $i$ and $i'$. We suppose that all of the units $i$ have a positive inclusion probability $\pi_i > 0$ and we denote by $d_i = 1/\pi_i$ their sampling weights. Two standard sampling designs are considered in the present paper : **Simple Random Sampling Without Replacement (SRSWOR)** of size $n_F$, and Poisson design with inclusion probabilities $\pi_i$, $i \in U_F$. For **SRSWOR**, $p_F$ assigns an equal probability to all without replacement samples of size $n_F$ and zero otherwise. The sampling weights are equal to $d_i = N_F/n_F$ for all $i$ in $U_F$. For Poisson sampling, the variables $I_i$'s are independent and distributed as Bernoulli random variables with parameter $\pi_i$.

The sample $s_F$ in $U_F$ leads to a sample $s_T$ in $U_T$, which is made of the units in $U_T$ linked to at least one unit in $s_F$. However, the sampling design $p_T(\cdot)$ which governs the selection of $s_T$, as well as the associated first-order inclusion probabilities, may be difficult to derive (Deville and Lavallée, 2006). Fortunately, as we will see in the next paragraph, for the **GWSM** estimators only $p_F(\cdot)$ and the associated inclusion probabilities are needed for the estimation of $t_y$ and $p_T(\cdot)$ will not be used.

Consider, for all $i \in U_F$ and $k \in U_T$, some non negative link weight $\theta_{ik}$ associated to the link $l_{ik}$ between $U_F$ and $U_T$, such that $\theta_{ik}$ is positive when $l_{ik} = 1$ and $\theta_{ik} = 0$ otherwise. We define the standardized link weights $\tilde{\theta}_{ik} = \theta_{ik} / \sum_{i' \in U_F} \theta_{i'k}$ which satisfy the constraint $\sum_{i \in U_F} \tilde{\theta}_{ik} = 1$. To compute the standardized link weights $\tilde{\theta}_{ik}$ for a given $k$ in $U_T$, one needs to know $\sum_{i \in U_F} \theta_{ik}$, which implies that the units $i$ in $U_F$ linked with $k$ must be known. We can take as an example $\theta_{ik} = l_{ik}$, and in this case, standardization implies that the number of units $i$ in $U_F$ linked with $k$ is known. More general weights can also be considered ; see Deville and Lavallée (2006) and Haziza and Beaumont (2017). The total $t_y$ can then be written as
the total on $U_F$ of the variable $\sum_{k \in U_T} \tilde{\theta}_{ik} y_k, i \in U_F$, as follows :

$$t_y = \sum_{k \in U_T} y_k = \sum_{k \in U_T} \left( \sum_{i \in U_F} \tilde{\theta}_{ik} \right) y_k = \sum_{i \in U_F} \left( \sum_{k \in U_T} \tilde{\theta}_{ik} y_k \right).$$

The estimation of $t_y$ can be obtained by considering standard estimators based on the sample $s_F$ selected from $U_F$. The Horvitz-Thompson (HT) estimator of $t_y$ 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) \qquad (2.1)$$

and it estimates unbiasedly the total $t_y$, provided that the link weights $\tilde{\theta}_{ik}, i \in U_F$, are standardized. This estimator is the Generalized Weight Share Method (GWSM) estimator, and it was studied by Deville and Lavallée (2006). It can also be written as follows :

$$\hat{t}_{y1} = \sum_{k \in U_T} \hat{t}_{\tilde{\theta}_k} y_k$$

where $\hat{t}_{\tilde{\theta}_k} = \sum_{i \in s_F} d_i \tilde{\theta}_{ik}$ is the HT estimator of the total $t_{\tilde{\theta}_k} = \sum_{i \in U_F} \tilde{\theta}_{ik} = 1$, for all $k$ in $U_T$. To calculate $\hat{t}_{y1}$, we only need to standardize the link weights that correspond to the sampled units $k$ in $U_T$.

Let us denote

$$\Delta_{ii'} = \frac{\pi_{ii'} - \pi_i \pi_{i'}}{\pi_i \pi_{i'}}, \quad i, i' \in U_F.$$

The variance of the GWSM estimator $\hat{t}_{y1}$ is given 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'}.$$

Interestingly, this variance can be rewritten as :

$$\mathrm{Var}(\hat{t}_{y1}) = \sum_{k \in U_T} \sum_{k' \in U_T} y_k y_{k'} \mathrm{Cov}(\hat{t}_{\tilde{\theta}_k}, \hat{t}_{\tilde{\theta}_{k'}}), \qquad (2.2)$$

where

$$\mathrm{Cov}(\hat{t}_{\tilde{\theta}_k}, \hat{t}_{\tilde{\theta}_{k'}}) = \sum_{i \in U_F} \sum_{i' \in U_F} \Delta_{ii'} \tilde{\theta}_{ik} \tilde{\theta}_{i'k'}.$$

The variance expression (2.2) is similar to the variance in the case of a direct sampling design on $U_T$. However, while for direct sampling, the covariance involved in the double sum is between the sample membership indicators weighted by the sampling weights, for indirect
sampling, the covariance is between the **HT** estimators of the link weights totals.

Deville and Lavallée (2006) were interested in finding the optimal weights $\tilde{\theta}_{ik}^{opt}$ that minimize the variance $\text{Var}(\hat{t}_{y1})$ for any survey variable $y$, namely

$$ \text{Var}(\hat{t}_{y1}^{opt}) \leq \text{Var}(\hat{t}_{y1}), \text{ for all } y, \qquad (2.3) $$

where

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

is the **GWSM** estimator obtained by using the optimal link weights $\tilde{\theta}_{ik}^{opt}$, $i \in U_F$, $k \in U_T$. Such a criterion is called *strong optimality* criterion and the optimal weights derived in this way, if they exist, should not depend on the survey variable $y$. However, the optimal link weights satisfying (2.3) might not exist, even for particular sampling designs such as Poisson or SRSWOR sampling designs (see Deville and Lavallée, 2006). Therefore, Deville and Lavallée (2006) suggested the *weak optimality* criterion which consists in finding the weak-optimal weights $\tilde{\theta}_{ik}^{wopt}$ which minimize the variance $\text{Var}(\hat{t}_{y1})$ for the particular variables $y$ such that $y_k = 1$ for a unit $k \in U_T$ and $y_{k'} = 0$ for $k' \neq k \in U_T$. We obtain for these particular variables :

$$ \hat{t}_{y1} = \sum_{i \in s_F} d_i \tilde{\theta}_{ik} = \hat{t}_{\tilde{\theta}_k}. $$

Thus, the weak-optimal weights $\tilde{\theta}_{ik}^{wopt}$ are the weights that minimize the variance of $\hat{t}_{\tilde{\theta}_k}$:

$$ (\tilde{\theta}_{ik}^{wopt})_{i \in U_F} = \arg \min_{\tilde{\theta}_{ik}, i' \in U_F} \text{Var}(\hat{t}_{\tilde{\theta}_k}), \quad \text{for all } k \in U_T. \qquad (2.4) $$

Deville and Lavallée (2006) derived weak-optimal link weights for Poisson and SRSWOR, and noticed that the weak optimality is a necessary condition for strong optimality.

### 2.2.2 MtO links and optimal weight links

The strong minimization problem previously mentioned becomes easier to handle when the links between $U_F$ and $U_T$ are Many-to-One. With MtO links, every unit from $U_F$ is linked to only one unit from $U_T$, and we can order the units in $U_F$ with respect to their common linked unit in $U_T$. For a given unit $k \in U_T$, let us denote by $U_{Fk}$, with size $N_{Fk} = \sum_{i \in U_F} l_{ik}$, the set of units $i$ in $U_F$ that are linked to $k$. In what follows, we consider that units in $U_F$ are ordered according to these subpopulations. Let $\Delta = (\Delta_{ii'})_{i,i' \in U_F}$ be the matrix of size $N_F \times N_F$. Thanks to this ordering, we can consider the submatrix $\Delta_{kk'} = (\Delta_{ii'})_{i \in U_{Fk}, i' \in U_{Fk'}}$ of $\Delta$ corresponding to elements in positions $i$ and $i'$ such that $i$ (resp. $i'$) is linked to $k$ (resp.
$k'$), for all $k$ (resp. $k'$) in $U_T$. For simplicity, we denote $\Delta_k$ the $\Delta_{kk}$ square submatrix with size $N_{Fk}$. With MtO links, the submatrices $\Delta_{kk'}$, $k, k' \in U_T$, form a partition of $\Delta$, namely

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

Let $1_k$ be the $N_{Fk}$-dimensional vector of ones. For MtO links, a sampling design is said to satisfy the $\Delta$-property if, for all $k \in U_T$, $\Delta_k$ is invertible and, for $k \neq k' \in U_T$, we have

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

The $\Delta$-property holds for Poisson sampling, SRSWOR and stratified SRSWOR under conditions detailed below.

For Poisson sampling from $U_F$ with inclusion probabilities $\pi_i$, $i \in U_F$, $\Delta_k$ is diagonal with positive terms, thus invertible, and $c_{kk'} = 0$ for all $k \neq k' \in U_T$. For SRSWOR of size $n_F$ from $U_F$, $\Delta_k$ is invertible as soon as $N_T > 1$, as can be seen in Deville and Lavallée (2006), page 174. If we denote $f = n_F/N_F$, we have

$$ c_{kk'} = -\frac{1-f}{f} \frac{1}{N_F-1}, $$

which does not depend on $i$ and $i'$, for all $k \neq k' \in U_T$. 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$. The submatrix $\Delta_{k,k' \neq k}$ cannot generally be written as $c_{kk'} 1_k 1_{k'}^t$, especially if $k$ is linked to one unit $i$ from stratum $h$ and one unit $i'$ from stratum $h' \neq h$. The $\Delta$-property holds if we assume that, for all $k$ in $U_T$, all units $i$ linked to $k$ belong to the same stratum $h$ and if, for each stratum $h$, there are at least two units of $U_T$ linked to $h$. Then $\Delta_k$ is invertible, and we have

$$ c_{kk'} = \begin{cases} -\frac{1-f_h}{f_h} \frac{1}{N_h-1}, & \text{if all units in } U_F \text{ linked to } k \text{ and } k' \text{ are in the same stratum } h \\ 0 & \text{otherwise.} \end{cases} $$

Thus, $c_{kk'}$ does not depend on $i$ and $i'$, for all $k \neq k' \in U_T$.

The first part of Proposition 2.2.1 gives an expression of the variance of the GWSM estimator while the second part gives the optimal link weights, in the MtO case, for sampling designs that satisfy the $\Delta$-property.

**Proposition 2.2.1.** If the links are MtO and the sampling design satisfies the $\Delta$-property,
then :

$$ Var(\hat{t}_{y1}) = \sum_{k \in U_T} y_k^2 Var(\hat{t}_{\tilde{\theta}_k}) + \sum_{k \in U_T} \sum_{k' \neq k \in U_T} y_k y_{k'} c_{kk'} \quad (2.6) $$

where $\hat{t}_{\tilde{\theta}_k} = \sum_{i \in s_F} d_i \tilde{\theta}_{ik}$ is the HT estimator of the total $t_{\tilde{\theta}_k} = 1$ for all $k \in U_T$.

Moreover, there is a unique set of optimal link weights verifying the strong-optimality criterion given in (2.3), given by

$$ (\tilde{\theta}_{ik}^{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. $$

The second term in the right-hand term of (2.6), does not depend on the link weights. As a consequence, minimizing the variance of $\hat{t}_{y1}$, regardless of the variable $y$ is, is equivalent to minimizing $Var(\hat{t}_{\tilde{\theta}_k})$ for all $k$ in $U_T$. Using the terminology introduced by Deville and Lavallée (2006) (see also relations (2.3) and (2.4)), this result means that, for MtO links and sampling designs that satisfy the $\Delta$-property, weak-optimality of the link weights is equivalent to strong-optimality and the link weights will be simply called optimal in the following.

For Poisson sampling, we have $\Delta_{ii} = (1 - \pi_i)/\pi_i$ and the optimal link weights are equal to :

$$ \tilde{\theta}_{ik}^{\text{opt}} = \frac{l_{ik}/\Delta_{ii}}{\sum_{i' \in U_F} l_{i'k}/\Delta_{i'i'}} \quad (2.7) $$

for all $i \in U_F$ and $k \in U_T$. For SRSWOR sampling, the optimal link weights are equal to $\tilde{\theta}_{ik}^{opt} = l_{ik}/\sum_{i' \in U_F} l_{i'k}$, for all $i \in U_F$ and $k \in U_T$. Details on the derivation of these optimal link weights can be found in Deville and Lavallée (2006). For stratified SRSWOR with the assumptions previously mentioned, it is easy to prove that the optimal link weights are the same as those for SRSWOR by following the proof in Deville and Lavallée (2006).

Consider now the GWSM estimator $\hat{t}_{y1}$ given in (2.1) and computed with some standardized link weights $\tilde{\theta}_{ik}$, $i \in U_F$, $k \in U_T$. Consider also the optimal GWSM estimator denoted by $\hat{t}_{y1}^{opt}$, and computed with the optimal link weights $\tilde{\theta}_{ik}^{opt}$. For MtO links and designs that satisfy the $\Delta$-property, it is possible to derive a new formula for the loss of efficiency between $\hat{t}_{y1}$ and $\hat{t}_{y1}^{opt}$. This loss can be expressed as a simple function of the variances of the HT estimator of $t_{\tilde{\theta}_k - \tilde{\theta}_k^{opt}} = \sum_{i \in U_F} (\tilde{\theta}_{ik} - \tilde{\theta}_{ik}^{opt}) = 0$, $k$ in $U_T$.

**Proposition 2.2.2.** If the links are MtO and the sampling design satisfies the $\Delta$-property, 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 :

As mentioned before, the matrices $\Delta$ derived with Poisson and SRSWOR designs satisfy
the Δ-property. Thus, we can compute easily the loss of efficiency between $\hat{t}_{y1}$ and $\hat{t}_{y1}^{opt}$ as formulated in the following corollary.

**Corollary 2.2.1.** If the links are MtO and if the sampling designs are Poisson or SRSWOR, then the loss in efficiency between $\hat{t}_{y1}$ and $\hat{t}_{y1}^{opt}$ has the expression given in Proposition 2.2.2. For Poisson sampling, we have:

$$Var(\hat{t}_{y1}) - Var(\hat{t}_{y1}^{opt}) = \sum_{k \in U_T} y_k^2 \sum_{i \in U_F} \frac{1 - \pi_i}{\pi_i} (\tilde{\theta}_{ik} - \tilde{\theta}_{ik}^{opt})^2.$$

For SRSWOR, we have:

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

$$\text{where } c = N_F^2 \left( \frac{1}{N_F} - \frac{1}{N_F} \right) \frac{1}{N_F - 1}.$$

The previous expression for SRSWOR can be easily adapted to stratified SRSWOR under the assumptions mentioned previously.

## 2.2.3 Comparison of direct and optimal indirect MtO sampling designs

It is not possible to compare theoretically the variances of simple indirect and direct estimators for general sampling designs even if we restrict ourselves to MtO links.

However, for Poisson sampling, we prove in Proposition 2.2.3 that the variance of the direct HT estimator is always smaller than the variance of the simple GWSM estimator when using optimal link weights. Let us consider a Poisson sampling design in $U_F$ with first-order inclusion probabilities $\pi_i$, $i \in U_F$. Because of the independence between the inclusion indicators, this sampling design induces a Poisson sampling design on $U_T$. Accounting for the MtO links between $U_F$ and $U_T$, we can calculate the probability of inclusion of every unit $k \in U_T$. We have that

$$\pi_k = P(\text{at least one } i \in U_F \text{ linked to } k \text{ is sampled in } s_F) = 1 - \prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}.$$

**Proposition 2.2.3.** Let us consider a sample $s_F$ drawn in a population $U_F$, using a Poisson sampling design with inclusion probabilities $0 < \pi_i < 1$. Let $U_T$ be another population associated to $U_F$ through MtO links $l_{ik}$. The sample $s_T$ deduced from $s_F$ using the MtO links between $U_F$ and $U_T$ can be considered as drawn in $U_T$, using Poisson sampling with inclusion probabilities $\pi_k = 1 - \prod_{i \in U_F} (1 - \pi_i)^{l_{ik}}$. The variance of the direct HT estimator, $\hat{t}_y = \sum_{k \in s_T} y_k / \pi_k$,
is smaller than the variance of the GWSM estimator $\hat{t}_{y1}^{opt} = \sum_{i \in s_F} \sum_{k \in U_T} \tilde{\theta}_{ik}^{opt} y_k / \pi_i$ calculated
with optimal link weights given in Equation (2.7).

As already stated, for Poisson sampling and MtO links, the optimal link weights lead to
the smallest possible variance of the simple GWSM estimator for any variable of interest.
So, under the assumptions of Proposition 2.2.3, the simple indirect estimator is always less
precise than the direct estimator.

The MtO case is interesting because it is possible, at least for Poisson sampling, to com-
pare the direct and the indirect sampling designs. Moreover, for several standard sampling
designs, it is possible to define optimal link weights and to calculate the exact loss of preci-
sion when using non-optimal link weights. However, when the number of units in the frame
population linked with a unit in the target population is large, all of the links might be not
observable, and an MtO indirect sampling could be very costly or even unfeasible. This pro-
blem arises for the La Poste survey in which the number of addresses per postman round is
500, on average, and where it is not possible to enumerate all addresses before the departure
of the postman in the morning. One solution is to use a double indirect sampling design as
detailed in the next section.

2.3 Double indirect sampling

2.3.1 Double GWSM

Double indirect sampling or indirect sampling in two steps (see Deville and Lavallée, 2006) consists of introducing an intermediate population $U_M$ in between the frame and the target populations, and using the same principle as for simple indirect sampling. There could be various reasons for introducing such a population. One reason could be that the target population $U_T$ units are only reachable through $U_M$. Another reason could be to simplify derivations. For example, Deville and Lavallée (2006) introduces an artificial intermediate population to simplify the search of an optimal standardized link matrix. In the La Poste survey, the objective is rather to decrease the number of links to observe. At La Poste, the intermediate population is a population of mail sorting boxes. Every morning, postmen sort the letters into boxes and deliver the letters from their allocated boxes (see Figure 3.2). The population of boxes is used as an intermediate population $U_M$ to link the addresses from the frame population $U_F$ to the postman rounds from the target population $U_T$.

Let $N_M$ be the size of the intermediate population $U_M$. Let $l_{ij}^{FM}$ be the link between $i \in U_F$ and $j \in U_M$ and let $l_{jk}^{MT}$ be the link between $j \in U_M$ and $k \in U_T$. A unit $i$ from the
<figure><img src="image_37.png" /><figcaption>FIGURE 2.2 – The frame, intermediate and target populations at La Poste.</figcaption></figure>

frame population $U_F$ could be linked to a unit $k$ from the target population $U_T$ by means of the unit $j$ from the intermediate population. The three populations $U_F$, $U_M$ and $U_T$ could be linked in various ways, as in simple indirect sampling with MtO, OtM and MtM links (see Figure 2.1).

As in Section 2.2, we consider non negative link weights $\theta_{ik}$ associated with the links $l_{ik}$ between $U_F$ and $U_T$ such that $\theta_{ik}$ is positive when $l_{ik} = 1$ and $\theta_{ik} = 0$ otherwise. We consider also the non negative weights $\theta_{ij}^{FM}$ associated with the links $l_{ij}^{FM}$ between $U_F$ and $U_M$ such that $\theta_{ij}^{FM}$ is positive when $l_{ij}^{FM} = 1$ and $\theta_{ij}^{FM} = 0$ otherwise. Finally, $\theta_{jk}^{MT}$ are non negative weights associated with the links between $U_M$ and $U_T$ and are defined in a similar way. With this double indirect sampling design, it can be seen that the links between units $i$ from $U_F$ and $k$ from $U_T$ are weighted by $\sum_{j \in U_M} \theta_{ij}^{FM} \theta_{jk}^{MT}$, and this link weight could be different from the link weight $\theta_{ik}$ used in simple indirect sampling.

Let $\tilde{\theta}_{ik}$ be the standardized link weight used in simple indirect sampling from $U_F$ to $U_T$, namely $\sum_{i \in U_F} \tilde{\theta}_{ik} = 1$ for all $k \in U_T$. We denote by $\tilde{\theta}_{ij}^{FM}$, respectively $\tilde{\theta}_{jk}^{MT}$, the link weights between $U_F$ and $U_M$, respectively $U_M$ and $U_T$, such that the link weights used in double indirect sampling between the frame population $U_F$ and the target population $U_T$ are normalized, namely $\sum_{i \in U_F} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} = 1$ for all $k \in U_T$. Note that this standardization does not require us to standardize each set of links $\tilde{\theta}_{ij}^{FM}$ and $\tilde{\theta}_{jk}^{MT}$. However, to obtain these standardized link weights, for a given $k$ in $U_T$, one needs to know the sum of the link weights of units from $U_F$ linked to $k$ passing by the intermediate population $U_M$, namely $\sum_{i \in U_F} \sum_{j \in U_M} \theta_{ij}^{FM} \theta_{jk}^{MT}$ should be known. The finite population total $t_y$ of $y$ can then be written as a total on the frame population $U_F$ as follows :

$$t_y = \sum_{k \in U_T} y_k = \sum_{k \in U_T} \left( \sum_{i \in U_F} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} \right) y_k = \sum_{i \in U_F} \left( \sum_{k \in U_T} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} y_k \right).$$
An estimator of $t_y$ can be derived easily by using the unbiased **HT** estimator as follows :

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

We call $\hat{t}_{y2}$ the double **GWSM** estimator, and its variance is given by :

$$ \mathrm{Var}(\hat{t}_{y2}) = \sum_{i,i' \in U_F} \Delta_{ii'} \sum_{k \in U_T} \sum_{j \in U_M} \hat{\theta}_{ij}^{FM} \hat{\theta}_{jk}^{MT} y_k \sum_{k' \in U_T} \sum_{j' \in U_M} \hat{\theta}_{i'j'}^{FM} \hat{\theta}_{j'k'}^{MT} y_{k'}. $$

### 2.3.2 MtO links

In this subsection, we focus on MtO links between the frame and the target populations. Comparing expressions (2.1) and (2.8), we can deduce that the double **GWSM** estimator $\hat{t}_{y2}$ can be viewed as a simple **GWSM** estimator with link weights $\hat{\theta}_{ik} = \sum_{j \in U_M} \hat{\theta}_{ij}^{FM} \hat{\theta}_{jk}^{MT}$.

Assuming that the sampling design verifies the **Δ**-property, we can apply Proposition 2.2.2 to determine the differences of the variances between double and simple **GWSM**. Thus, Proposition 2.2.2 shows that, for any variable of interest $y$, the optimal simple **GWSM** estimator is always better than the double **GWSM** estimator. The loss of efficiency of the double **GWSM** estimator with respect to the optimal simple **GWSM** estimator depends on the variances of the **HT** estimators $\sum_{i \in s_F} d_i(\sum_{j \in U_M} \hat{\theta}_{ij}^{FM} \hat{\theta}_{jk}^{MT} - \hat{\theta}_{ik}^{opt})$ for $k \in U_T$. This loss depends on the configuration of the link weights used in the double indirect sampling. If the double indirect sampling weights, $\sum_{j \in U_M} \hat{\theta}_{ij}^{FM} \hat{\theta}_{jk}^{MT}$, are close to the optimal simple indirect sampling weights, $\hat{\theta}_{ik}^{opt}$, for all $i \in U_F$, then the use of a double **GWSM** estimator will cause a small loss in precision. Otherwise, the loss could be substantial.

In the following subsection, we describe a particular double indirect sampling design that requires fewer links than its simple indirect sampling counterpart while maintaining the same precision.

### 2.3.3 MtO-MtO links and double standardization

Consider the MtO-MtO case where the links between the frame and the intermediate population are MtO, and the links between the intermediate and the target population are also MtO. This case implies that the links between $U_F$ and $U_T$ are MtO. In this situation, the double **GWSM** estimator has a simple expression, since a unit $i$ from the frame population is linked to a single unit $j$ from the intermediate population, which is itself linked to a single unit $k$ from the target population. Thus, for a given sampled unit $i \in U_F$, the sums over the intermediate and the target populations in the double **GWSM** estimator $\hat{t}_{y2}$ given by (2.8),
<figure><img src="image_38.png" /><figcaption>FIGURE 2.3 – Number of observed links for each standardization</figcaption></figure>

contain only one non-zero element equal to $\tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} y_k$. To compute $\hat{t}_{y2}$, we need to compute only the standardized link weight $\tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}$ that corresponds to the unique unit $k \in U_T$ linked to the sampled unit $i \in U_F$ through the unique $j \in U_M$.

In this MtO-MtO setup, the choice of the standardization method arises and has an impact on the number of links to observe. One can consider link weights $\tilde{\theta}_{ij}^{FM}$ and $\tilde{\theta}_{jk}^{MT}$ that are either both standardized ($\sum_{i' \in U_F} \tilde{\theta}_{i'j}^{FM} = 1$ and $\sum_{j' \in U_M} \tilde{\theta}_{j'k}^{MT} = 1$), or globally standardized ($\sum_{i' \in U_F} \sum_{j' \in U_M} \tilde{\theta}_{i'j}^{FM} \tilde{\theta}_{j'k}^{MT} = 1$). Note that if both link weights are standardized (double standardization), then they are also globally standardized, but the converse is not true. As detailed below, using MtO-MtO with double standardization could allow for a reduction in the number of links to observe compared to the MtO-MtO with global standardization, or even compared to the simple MtO GWSM.

With double standardization, we can derive $\tilde{\theta}_{ij}^{FM}$ and $\tilde{\theta}_{jk}^{MT}$ separately. To compute $\tilde{\theta}_{ij}^{FM}$, the links between $U_F$ and the units $j \in U_M$ indirectly sampled through $i \in s_F$ need to be observed, which gives the total number of links to observe equal to $N_{Fj} = \sum_{i' \in U_F} l_{i'j}^{FM}$. Similarly, to compute $\tilde{\theta}_{jk}^{MT}$, the links between $U_M$ and the indirectly sampled $k \in U_T$ must be observed which leads to $N_{Mk} = \sum_{j' \in U_M} l_{j'k}^{MT}$ links to observe. Thus, for MtO-MtO GWSM with double standardization and for a given unit $i \in s_F$, the total number of links to observe is equal to $N_{Fj} + N_{Mk}$. The right plot of Fig. 2.3 shows the links to observe (black links) on a small example for MtO-MtO with double standardization.

For MtO-MtO with global standardization, we need to observe the links between $U_F$ and the units $j \in U_M$ linked to the indirectly sampled units $k \in s_T$, namely we need to know $\sum_{i' \in U_F} \sum_{j' \in U_M} l_{i'j}^{FM} l_{j'k}^{MT} = \sum_{j' \in U_M} N_{Fj'} l_{j'k}^{MT}$. We also need the number of links between $U_M$ and the indirectly sampled $k \in U_T$, namely we need to know $N_{Mk} = \sum_{j' \in U_M} l_{j'k}^{MT}$. Thus, for MtO-MtO links with global standardization, we need to observe $\sum_{j' \in U_M} N_{Fj'} l_{j'k}^{MT} + N_{Mk}$
links. The middle plot of Fig. 2.3 shows the links to observe (black links), for a given unit $i \in s_F$, for MtO-MtO with global standardization.

Consider now simple MtO indirect sampling. A unit $i$ from the frame population can be linked to a single unit $k$ from the target population, and the sum over the target population in the simple GWSM estimator $\hat{t}_{y1}$ given by (2.1) contains only a non-zero element equal to $\tilde{\theta}_{ik}y_k$. To compute $\hat{t}_{y1}$, for each unit $i \in s_F$, we need to compute the standardized link weight $\tilde{\theta}_{ik}$ that corresponds to the unique unit $k \in U_T$ linked to $i \in s_F$. This circumstance implies that we need to observe the links between the indirectly sampled unit $k$ of the target population and the frame population; namely we need to observe a number of links equal to $N_{Fk} = \sum_{i' \in U_F} l_{i'k}$. The left plot of Fig. 2.3 shows the links to observe (black links) for a given unit $i \in s_F$ in MtO indirect sampling.

It is possible to compare the number of links between MtO-MtO with double standardization with simple MtO if we assume that $N_{Fk} = N_{Fj}N_{Mk}$, $N_{Fj} > 2$ and $N_{Mk} > 2$. We then have

$$N_{Fj} + N_{Mk} < N_{Fj}N_{Mk} = N_{Fk}, \qquad (2.9)$$

and the double GWSM with double standardization always requires fewer links to observe than the simple GWSM. Moreover, the smallest number of links to observe with the double GWSM with double standardization is achieved when $N_{Fj} = N_{Mk} = N_{Fk}^{1/2}$, which is the most favorable situation for the double GWSM.

The assumption $N_{Fk} = N_{Fj}N_{Mk}$ is true if the numbers of links $N_{Fj}$ are equal for all $j$ linked to the same $k$ in $U_T$. Indeed, if we let $C_k$ denote a positive constant, that does not depend on $j$, such that $N_{Fj} = C_k$ for all units $j \in U_M$ linked to $k$ in $U_T$, then,

$$N_{Fk} = \sum_{j \in U_M} l_{jk}^{FM} N_{Fj} = C_k \sum_{j \in U_M} l_{jk}^{MT} = C_k N_{Mk} = N_{Fj} N_{Mk}. \qquad (2.10)$$

This remark is interesting as it stands in that when the number of links between $U_F$ and $U_M$ is the same for each unit in $U_T$, there will be fewer links to observe when using MtO-MtO with double standardization, regardless of what are the links between $U_M$ and $U_T$ (see Setups 1 and 3 in Section 2.4).

Furthermore, if the link weights are the link indicators, namely $\theta_{ij}^{FM} = l_{ij}^{FM}$, $\theta_{jk}^{MT} = l_{jk}^{MT}$, $\theta_{ik} = l_{ik}$, and $N_{Fk} = N_{Fj}N_{Mk}$, then the double GWSM and the simple GWSM estimators are equal. Indeed,

$$\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} = \sum_{j \in U_M} \frac{l_{ij}^{FM} l_{jk}^{MT}}{N_{Fj} N_{Mk}} = \frac{l_{ik}}{N_{Fk}} = \tilde{\theta}_{ik}. \qquad (2.11)$$

In such a situation, the double GWSM with double standardization and the simple GWSM have the same precision, but the double GWSM ensures a gain in terms of the smaller number
of links to observe. However, if the condition $N_{Fk} = N_{Fj}N_{Mk}$ is not fulfilled, (see Setups 2 and 4 in Section 2.4 for further details), the double GWSM could be less efficient than the simple GWSM estimator.

For the double GWSM with global standardization, the number of links to observe is $\sum_{j' \in U_M} l_{j'k}^{MT} N_{Fj'} + N_{Mk} = N_{Fk} + N_{Mk}$, which is greater than the number of links to observe for both, the simple GWSM and the double GWSM with double standardization.

In Fig. 2.3, $N_{Fk} = 10$, $N_{Fj} = 5$ and $N_{Mk} = 2$. The number of links to observe is 10 for the simple GWSM (left plot of the figure), 12 for the double GWSM with global standardization (middle plot) and 7 for the double GWSM with double standardization (right plot). In the La Poste situation, where the double GWSM with double standardization is used, the gain is much larger because on average $N_{Fk} = 500$, $N_{Fj} = 10$ and $N_{Mk} = 50$. Thus, we have $N_{Fj} + N_{Mk} = 60$ links to observe for the double GWSM with double standardization while it is $N_{Fj}N_{Mk} = 500$ for the simple GWSM and $N_{Fk} + N_{Mk} = 550$ for the double GWSM with global standardization.

Double GWSM with global standardization is of poor interest in the context of La Poste. The loss of precision is not compensated by a gain in the number of links to observe. In what follows, we focus only on double GWSM with double standardization.

## 2.4 Simulation study

We have shown in Subsection 2.3.1 that, in the case of MtO links and sampling designs with the Δ-property, the loss of precision of the double GWSM compared to the optimal simple GWSM depends on the variance of the HT estimators $\sum_{i \in s_F} d_i(\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT} - \tilde{\theta}_{ik}^{opt})$ for all $k \in U_T$. In other words, the increase in the variance depends on the configuration of the double GWSM link weights $\sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \tilde{\theta}_{jk}^{MT}$, compared to the optimal simple link weights $\tilde{\theta}_{ik}^{opt}$. In this section, we conduct a Monte Carlo study to analyze the influence of the link weights on the efficiency of the double GWSM estimator.

### 2.4.1 Population and link setups

We have generated three populations, $U_F, U_M$ and $U_T$ as well as the links between them, to meet the framework assumed in Subsection 2.3.1 as well as at the La Poste situation. At La Poste, the target population $U_T$ is the population of rounds, the frame population $U_F$ is the population of addresses, and the intermediate population $U_M$ is the population of boxes. All links are MtO, which means that an address can only be in one box and a box in only one round. There are on average 50 boxes per round and 500 addresses per round. The three
populations $U_F, U_M$ and $U_T$ were generated in the following way. The target population of rounds $U_T$ of size $N_T = 6958$ was obtained from La Poste's historical data. Then, based on $U_T$, we generated the intermediate population of boxes $U_M$ of size $N_M = 347900 = 50 \times 6958$, and the frame population of addresses $U_F$, of size $U_F = 3479000 = 500 \times 6958$, together with the links.

We are interested in estimating the total $t_y = \sum_{k \in U_T} y_k$ of a study variable $y$ which is a particular measure of the postal traffic obtained from La Poste's historical data. For confidentiality reasons, the study variable $y$ was transformed.

We considered different setups of MtO links between $U_F$ and $U_M$ and between $U_M$ and $U_T$. For ease of comparison, we kept unchanged the links between $U_F$ and $U_T$ in all scenarios. In this double indirect sampling design, additional care must be taken to ensure that, if unit $i \in U_F$ is linked to $k \in U_T$ in simple indirect sampling case, then it must be linked to $k \in U_T$ for double indirect sampling also. There is a link between $U_F$ (resp. $U_F, U_M$) and $U_T$ (resp. $U_M, U_T$) when an address (resp. an address, a box) is part of a round (resp. a box, a round). The first set of links was generated between $U_F$ and $U_T$, and is such that the number of addresses $N_{Fk}$ in a round $k$ of $U_T$ is equal to 500, for all rounds $k \in U_T$. Thus, $N_{Fk} = \sum_{j \in U_M} l_{jk}^{MT} N_{Fj} = 500$, where $N_{Fj}$ is the number of addresses in the box $j$. This configuration means that all units from the target population have the same number of links with the frame population. As detailed in Subsection 2.4.2, the advantage of this simplification is that, under some supplementary assumptions on the populations and on the inclusion probabilities that are true in our setting, the simple GWSM estimator has the same expression as a direct HT estimator. To study the effect of the links structure on the efficiency of the double GWSM estimator, we created four different setups of links with the intermediate population. The links between $U_F$ and $U_M$ (resp. $U_M$ and $U_T$) were generated either uniformly or not. The links between $U_F$ and the boxes $j$ of $U_M$ are called “uniform” when there is the same number of addresses, $N_{Fj}$, in all boxes $j$ of $U_M$ that are part of the same round $k$ of $U_T$. They are called “non-uniform” when we generate one address in all boxes of $U_M$ that are part of the same round of $U_T$, except for one box $j_0$ which contains the remaining addresses. This last situation is a type of extreme unbalanced case for the number of links between $U_F$ and $U_M$. The links between $U_M$ and $U_T$ are “uniform” when there is the same number of boxes $N_{Mk}$ per round. The “non-uniform” case is generated by considering two boxes for 6286 rounds and 499 boxes for the remaining 672 rounds. Note that we cannot choose non-uniform links between $U_M$ and $U_T$ with one box or 500 boxes in a round together with non-uniform links between $U_F$ and $U_M$ similar to the ones proposed above. The reason is that we have set a constraint of 500 addresses per round. Under this constraint, having one box (resp. 500 boxes) in a round implies having 500 addresses (resp. one address) in each
box of the round, which corresponds to uniform links between the addresses and the boxes.

The four setups are detailed below (see also Figure 2.4 for graphical examples) :

- Setup 1 : Uniform/Uniform : the number of links between $U_F$ (resp. $U_M$) and $U_M$ (resp. $U_T$) is uniform. In this setup, the double GWSM and the simple GWSM estimators are equal, as proved in Subsection 2.3.3 (see (2.11)).

- Setup 2 : non-uniform/Uniform : the links between $U_F$ and $U_M$ are non-uniform while the links between $U_M$ and $U_T$ are uniform.

- Setup 3 : Uniform/non-uniform : the links between $U_F$ and $U_M$ are uniform while the links between $U_M$ and $U_T$ are non-uniform. In this setup, the $N_{Fj}$ must be rounded as illustrated on Figure 2.4 (bottom left panel) where, for the second unit in $U_T$, there are 6 addresses to divide between 4 boxes, and thus, 2 boxes contain 1 address each, while the other 2 boxes contain 2 addresses each. Ignoring the rounding of the $N_{Fj}$, the relation $N_{Fk} = N_{Fj}N_{Mk}$ holds (see also equation (2.10)), which allows for the equality of the simple and double GWSM estimators.

- Setup 4 : non-uniform/non-uniform : the links between $U_F$ and $U_M$, and between $U_M$ and $U_T$ are non-uniform.

## 2.4.2 Sampling designs and GWSM estimators

We consider two sampling designs that satisfy the $\Delta$-property : the SRSWOR of sizes $n = 500$ and $n = 1000$, as well as the Bernoulli design (which is a Poisson design with equal inclusion probabilities) with expected sample sizes equal to 500 and 1000.

Let $y_k, k \in U_T$, be a measure of the postal traffic in round $k$, as mentioned at the beginning of section 2.4.1. For each link setup, we compare the double GWSM estimator to the simple GWSM estimator, both computed on samples $s_F$ drawn from the frame population $U_F$ of addresses :

$$ \hat{t}_{y1} = \sum_{i \in s_F} \frac{1}{\pi_i} \sum_{k \in U_T} \tilde{\theta}_{ik} y_k, \quad \hat{t}_{y2} = \sum_{i \in s_F} \frac{1}{\pi_i} \sum_{j \in U_M} \tilde{\theta}_{ij}^{FM} \sum_{k \in U_T} \tilde{\theta}_{jk}^{MT} y_k. $$

For the simple GWSM estimator, we consider the optimal link weights, which are equal to $\tilde{\theta}_{ik} = l_{ik}/N_{Fk}$ for SRSWOR and Bernoulli samplings. For the double GWSM, we consider $\tilde{\theta}_{ij}^{FM} = l_{ij}^{FM}/N_{Fk}$ and $\tilde{\theta}_{jk}^{MT} = l_{jk}^{MT}/N_{Mk}$.

The above simulation setting facilitates the comparison between the double GWSM estimator and the simple GWSM estimator by ensuring that only the double GWSM varies, while the simple GWSM remains fixed. In fact, the setting makes the simple GWSM very close to the direct HT estimator in the 4 setups. We compare also the double GWSM estimator to the direct HT estimator as calculated from samples $s_T^*$ drawn directly from the
<figure><img src="image_39.png" /></figure>

**Setup 1**

Uniform number of adresses in a box

Uniform number of boxes in a round

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

**Setup 2**

Non uniform number of adresses in a box

Uniform number of boxes in a round

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

**Setup 3**

Uniform number of adresses in a box

Non uniform number of boxes in a round

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

**Setup 4**

Non uniform number of adresses in a box

Non uniform number of boxes in a round

<caption>FIGURE 2.4 – The four setups</caption>
target population $U_T$ of rounds using **SRSWOR** of size $n_T$:

$$ \hat{t}_{HT} = \sum_{k \in s_T^*} \frac{y_k}{\pi_k} $$

If the $\pi_i$ are small, the probability of drawing two or more units from $U_F$ that are linked to the same $k$ in $U_T$ is small and we can approximate $\pi_k$ by $N_{Fk}\pi_i$. In this situation, every unit from the target population is generally linked with at most one unit from the sample $s_F$, and thus, for all $k \in s_T$, we have $\sum_{i \in s_F} l_{ik} \approx 1$. This hypothesis is true at La Poste because given the large number of postal addresses in the frame population, we never sample several addresses from the same round. Thus, we have:

$$ \begin{aligned} \hat{t}_{y1} &= \sum_{i \in s_F} \frac{1}{\pi_i} \sum_{k \in U_T} \tilde{\theta}_{ik} y_k = \sum_{k \in U_T} y_k \sum_{i \in s_F} \frac{1}{\pi_i} \tilde{\theta}_{ik} \\ &\approx \sum_{k \in s_T} y_k \frac{N_{Fk} \sum_{i \in s_F} l_{ik}}{\pi_k N_{Fk}} \approx \sum_{k \in s_T} \frac{y_k}{\pi_k} \end{aligned} $$

which means that, in our setting, the simple **GWSM** estimator is approximately equivalent to the direct **HT** estimator for the sample $s_T$. Since this formula does not depend on the simulation setup, the simple **GWSM** does not vary between the setups.

### 2.4.3 Measures of efficiency and results

We have considered $R = 100 000$ samples according to the **SRSWOR** and Bernoulli sampling designs with sizes or expected sizes 500 and 1 000. We have computed the Monte Carlo relative bias of the simple **GWSM** and the double **GWSM** estimators:

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

where $\hat{t}_{y1}^{(r)}$ and $\hat{t}_{y2}^{(r)}$ are the simple and the double **GWSM** estimates respectively, computed for the $r$-th sample, $r = 1, \dots, R$. As a measure of efficiency, we have computed the Relative Root Mean Square Error (RRMSE) of $\hat{t}_{y2}$ with respect to $\hat{t}_{y1}$, and with respect to $\hat{t}_{HT}$:

$$ RRMSE_{MC}(\hat{t}_{y1}) = \sqrt{\frac{MSE_{MC}(\hat{t}_{y2})}{MSE_{MC}(\hat{t}_{y1})}} \quad \text{and} \quad RRMSE_{MC}(\hat{t}_{HT}) = \sqrt{\frac{MSE_{MC}(\hat{t}_{y2})}{MSE_{MC}(\hat{t}_{HT})}} $$

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, $$
<table><caption>TABLE 2.1 – Relative bias and efficiency of the DGWSM estimate under different links setups.</caption><thead><tr><th>Design</th><th>Setup</th><th>n</th><th>RB<sub>MC</sub>(t̂<sub>y2</sub>)</th><th>RRMSE<sub>MC</sub>(t̂<sub>y1</sub>)</th><th>RRMSE<sub>MC</sub>(t̂<sub>HT</sub>)</th></tr></thead><tbody><tr><td rowspan="4">SRSWOR<br>100 000<br>simulations</td><td>Setup 1</td><td>500<br>1000</td><td>0.15<br>0.09</td><td>100.00<br>100.00</td><td>103.96<br>109.78</td></tr><tr><td>Setup 2</td><td>500<br>1000</td><td>-0.53<br>0.44</td><td>337.89<br>334.51</td><td>347.14<br>365.85</td></tr><tr><td>Setup 3</td><td>500<br>1000</td><td>0.07<br>0.12</td><td>100.01<br>100.02</td><td>105.35<br>106.12</td></tr><tr><td>Setup 4</td><td>500<br>1000</td><td>2.67<br>0.33</td><td>1173.36<br>1093.90</td><td>1216.95<br>1198.97</td></tr></tbody></table>

and $MSE_{MC}(\hat{t}_{y1})$, and $MSE_{MC}(\hat{t}_{HT})$ are defined similarly.

Table 2.1 contains the simulation results for the SRSWOR design. Similar results were obtained for Bernoulli sampling and are not reported here. The results are also very comparable for both sampling sizes. As expected, both GWSM estimators have a low Monte Carlo relative bias in all setups.

Moreover, we observe small differences between the RRMSE of the double compared to the simple GWSM, and the double GWSM compared to the direct HT estimators in all setups. This result was expected since the simple GWSM expression corresponds to a direct HT expression (see details at the end of Subsection 2.4.2). Small differences arise because the samples $s_T^*$, drawn directly in the target population, differ from the samples $s_T$ that are obtained through the samples $s_F$ using the links between $U_F$ and $U_T$.

For Setups 1 and 3, as shown in Subsection 2.3.3, $\hat{t}_{y1} = \hat{t}_{y2}$. Thus, $RRMSE_{MC}(\hat{t}_{y1}) = 100\%$ and $MSE_{MC}(\hat{t}_{y1}) = MSE_{MC}(\hat{t}_{y2})$. The small loss of precision between the GWSM estimators in Setup 3 occurs because the relation $N_{Fk} = N_{Fj}N_{Mk}$ is not exactly satisfied due to rounding errors.

For Setups 2 and 4, the equation $N_{Fk} = N_{Fj}N_{Mk}$ does not hold at all. In Setup 2, if the sampled address $i$ is alone in the box $j$, then $N_{Fj}N_{Mk} = 1 * 50$ which is far from $N_{Fk} = 500$. If $i$ is in the box containing 451 addresses, then $N_{Fj}N_{Mk} = 451 * 50$ which is also far from $N_{Fk}$. In Setup 4, the difference between $N_{Fk}$ and $N_{Fj}N_{Mk}$ is even larger because we also let the $N_{Mk}$ vary. We note an important loss of precision of the double GWSM estimator compared to the simple GWSM estimator in Setups 2 and 4.

The precision of the double GWSM estimator depends on how close the values $N_{Fj}N_{Mk}$ and $N_{Fk}$ are, for every $j$ linked to the same $k$. As proved in equation (2.10), the uniform link structure between $U_F$ and $U_M$ implies that $N_{Fk} = N_{Fj}N_{Mk}$, regardless of the link structure between $U_M$ and $U_T$. This remark helps to explain the good results for Setups 1 and 3 and
the poor results for Setups 2 and 4.

It is also interesting to compare the number of links to observe for the two GWSM estimators in each setup. The simple GWSM estimator requires the observation of 500 links per sampled round in each setup. To compute the double GWSM estimator, there are 60 links to observe on average per sampled round in Setup 1, 457 in Setup 2, 276 in Setup 3 and 500 in Setup 4 (the averages are rounded values). For Setups 1 and 3, the equation $N_{Fk} = N_{Fj}N_{Mk}$ (almost) holds, and thus, there will always be a gain in the number of links to observe (see equation (2.9)). It can be noted that this gain is even larger if the $N_{Mk}$ are uniform. For Setup 2, the gain in the number of links is limited, while there is no gain in Setup 4.

The simulations illustrate that the link structure between the three populations has a large impact on the double GWSM estimator in terms of the precision, but also in terms of the number of links to observe. In the ideal situation of Setup 1, there is a clear advantage of using double indirect sampling for an MtO-MtO situation, while it is not at all recommended in situations like the one illustrated in Setup 4.

## 2.5 Application to the French post data

Before 2008, La Poste sampled directly the postmen rounds to estimate the monthly postal traffic. After a reorganization of the post offices in 2008, the population of rounds became incomplete and La Poste had to use indirect sampling through the frame population of addresses. Because of the large number of links to observe, a simple GWSM estimator was not possible. La Poste had to consider a double MtO-MtO indirect sampling design and a double GWSM estimator with double standardization. The use of double indirect sampling, compared to the previous direct sampling method, led to a precision loss of the estimators. The estimated standard deviations of the estimators were increased by a factor between 2 and 3. To complete the theoretical results of Sections 2.2 and 2.3, we propose, through simulations and in a setup similar to La Poste, to evaluate the loss of precision due to using double indirect sampling, and to check if the calculated loss is of the same order as the loss observed in reality.

In this application, we focus on simple GWSM and double GWSM, which are both computed on a sample of addresses, and a direct HT estimator computed on a sample of rounds. The samples are drawn according to SRSWOR designs. The sampling design at La Poste is more complex and involves a stratification based on a typology of the post offices. This stratification of the postal addresses ensures that a round cannot belong to two different strata of addresses, and that every stratum contains at least two rounds. Thus,
<table><caption>TABLE 2.2 – Relative bias and comparison of RMSE, in percentages, for the double (D) GWSM, the simple (S) GWSM and the direct estimates in a setup comparable to La Poste.</caption><thead><tr><th>Design</th><th>n</th><th>GWSM</th><th>RB<br>y = 1</th><th>RRMSE<br>y = 1<br>rel. to S</th><th>RB<br>y = traffic</th><th>RRMSE<br>y = traffic<br>rel. to S</th><th>RRMSE<br>y = traffic<br>rel. to direct</th></tr></thead><tbody><tr><td>SRSWOR</td><td>500</td><td>S</td><td>-0.01</td><td>100</td><td>0.02</td><td>100</td><td>137.96</td></tr><tr><td>100 000</td><td>500</td><td>D</td><td>-0.03</td><td>219.00</td><td>0.08</td><td>166.69</td><td>229.98</td></tr><tr><td>simulations</td><td>1000</td><td>S</td><td>0.00</td><td>100</td><td>-0.01</td><td>100</td><td>144.55</td></tr><tr><td></td><td>1000</td><td>D</td><td>0.00</td><td>220.89</td><td>-0.07</td><td>163.53</td><td>236.39</td></tr></tbody></table>

the Δ-property holds for this sampling design (see details in Subsection 2.2.2 on stratified SRSWOR). The estimators are also more complex and involve calibration and winzorisation. Considering such complex designs and estimators is beyond the scope of the present study, and we focus on the SRSWOR sampling design with direct HT, and with simple and double indirect GWSM estimators.

We do not only look at the loss caused by the use of a double GWSM compared to a simple GWSM. We also examine the loss caused by the use of double indirect sampling compared to a direct sampling. The objective is to capture the total loss in the precision observed at La Poste when changing their sampling design from direct to doubly indirect.

The setup for the simulations below is close to the La Poste setup in the sense that the number of addresses in a box, and the number of boxes in a round were generated using observed distributions from La Poste data. Compared to the four setups in Section 2.4, the number of addresses in a round was not fixed at 500, but computed using the number of addresses in a box and the number of boxes in a round. The number of addresses in a box varies from 1 to 29 with two modes at 1 and 13, with rare observations at approximately 40 and 120. The number of boxes in a round varies from 28 to 73 with two modes at approximately 35 and 70, and rare observations between 100 and 1000. The average number of addresses in a box in this setup is 14, the average number of boxes in a round is 60 and the average number of addresses in a round is 841. This setup is close to Setup 4 in Section 2.4, but it has less variability in the number of links, with the number of addresses between 30 and 70 and the number of boxes between 1 and 29, while the number of addresses and boxes varies between 2 and 499 in Setup 4.

We consider two study variables $y$, and we are interested in estimating their totals on the target population of rounds. The first study variable is equal to 1 for all units, which gives a total over the target population equal to $N_T$, while the second variable is a confidential measure of postal traffic obtained from La Poste data. For the selection of indirect sample, we consider simple random sampling without replacement with respective sizes of 500 and
1000 selected from the frame population of addresses. For the selection of direct sample, We also consider a simple random sampling without replacement with respective sizes 500 and 1000 selected in the target population of rounds.

As in Section 2.4, we compute, for $R = 100000$ simulations, the Monte Carlo relative bias (RB) as a percentage of the simple (S) and double (D) **GWSM** estimators together with their mean square errors. For both variables ($y = 1$ and $y = \text{traffic}$), we use the **RRMSE** in percentage, to compare the double to the simple **GWSM** (see the two **RRMSE** columns relative to S in Table 2.2). For $y = \text{traffic}$, we also compare the double **GWSM** to the direct estimator (see the **RRMSE** relative to direct in the last column of Table 2.2). We note that, for $y = 1$, the MSE of the direct estimator is zero since the estimator is calibrated on the size of the population. Thus, the **RRMSE** of the double **GWSM** compared to the direct is infinite and not reported. We notice that there is almost no difference between the results obtained for the two sample sizes. As expected, given that all estimators are unbiased, the relative biases of the simple and double **GWSM** estimators in Table 2.2 are small for both variables of interest. Moreover, we observe a loss in precision by using the double **GWSM** estimator instead of the simple **GWSM** estimator for both variables. This loss is less important than those observed in Setups 2 and 4, since the differences in the weights between the double **GWSM** and the simple **GWSM** are smaller here and have less variability than in Setups 2 and 4. In Table 2.1 (see the last two columns), for Setup 2 (resp. Setup 4), the standard deviations are multiplied by a factor between 3 and 4 (resp. 10 and 13). In Table 2.2, the standard deviations are multiplied by a factor between 1 and 2 (resp. 2 and 3) for the “traffic” (resp. 1) variable when comparing D and S. The additional loss of precision between the S and the direct **HT** estimators is not negligible, and gives a factor between 1 and 2 for the “traffic” variable (1.38 for the sample size $n_F=500$ and 1.45 for $n_F = 1000$). In total, the standard deviations increase by a factor between 2 and 3 (approximately 2.3) when changing from direct to double indirect sampling. Interestingly, the loss of precision that we observe is of the same order as the loss observed in practice at La Poste.

## 2.6 Conclusion

The **MtO** situation in indirect sampling allows us to obtain optimal link weights for some classical sampling designs such as simple random sampling without replacement, Poisson sampling and stratified **SRSWOR**. In this context, it is also possible to derive an exact expression for the loss of precision when the link weights are not optimal. This expression shows that the increase of variance of the **GWSM** estimator depends on how far the link weights estimators are from the optimal link weights estimators. When the number of links
to observe is large, it is possible to introduce a double indirect sampling design that allows us to reduce the number of links to observe when using a double standardization. As illustrated by our simulations, the double GWSM with double standardization proves to be especially interesting in some specific cases. It allows for a reduction in the number of observed links while maintaining the level of precision of a simple GWSM. However, it can be less useful in other cases, with a considerable loss of precision and not an important reduction of the number of links to observe. In the La Poste situation, there is a clear reduction in the number of addresses to observe per round, but at the cost of a large loss of precision. One perspective to improve on the precision of the estimators at La Poste is to keep the double indirect sampling design but use simple indirect GWSM and to predict the unobserved link weights as proposed for example by Xu and Lavallée (2009) and Falorsi et al. (2019). Indeed, double indirect sampling helps on saving costs and should be maintained. However, the use of the double GWSM with double standardization may lead to a significant loss of precision. Thus, a perspective that is currently under study at La Poste is to predict the number of unobserved links that are needed for the simple GWSM by using a model along with auxiliary information available at the level of the intermediate population of boxes.

## 2.7 Appendix

### Proof of Proposition 2.2.1

We compute first the variance of the simple GWSM assuming MtO links and a sampling design that satisfy the Δ-property. We follow Deville and Lavallée (2006) and use matrix notations. Let $\tilde{\Theta}_k = (\tilde{\theta}_{ik})_{i \in U_F}$ be the $N_F$-dimensional vector of standardized link weights. With MtO links, we can write $\tilde{\Theta}_k$ as follows :

$$ \tilde{\Theta}_k = \begin{pmatrix} 0 \\ \tilde{\theta}_k \\ 0 \end{pmatrix}, \quad k \in U_T, \qquad (2.12) $$

where $\tilde{\theta}_k = (\tilde{\theta}_{ik})_{i \in U_F}$ is the $N_{Fk}$-dimensional vector of positive weighted links with $N_{Fk} = \sum_{i \in U_F} l_{ik}$, the number of units $i$ in $U_F$ linked to $k$ in $U_T$. The variance of the GWSM estimator
$\hat{t}_{y1}$ given in (2.2) can be written as follows :

$$
\begin{align*}
\mathrm{Var}(\hat{t}_{y1}) &= \sum_{k \in U_T} \sum_{k' \in U_T} y_k y_{k'} \mathrm{Cov}(\hat{t}_{\tilde{\theta}_k}, \hat{t}_{\tilde{\theta}_{k'}}) \\
&= \sum_{k \in U_T} \sum_{k' \in U_T} y_k y_{k'} \tilde{\Theta}_k^t \Delta \tilde{\Theta}_{k'} \\
&= \sum_{k \in U_T} \sum_{k' \in U_T} y_k y_{k'} \tilde{\theta}_k^t \Delta_{kk'} \tilde{\theta}_{k'}, \tag{2.13}
\end{align*}
$$

where $\mathbf{\Delta} = \left( \frac{\pi_{ii'} - \pi_i \pi_{i'}}{\pi_i \pi_{i'}} \right)_{i,i' \in U_F}$ and $\mathbf{\Delta}_{kk'} = (\Delta_{ii'})_{i \text{ linked to } k, i' \text{ linked to } k'}$ is the submatrix of $\mathbf{\Delta}$ of size $N_{Fk} \times N_{Fk'}$ corresponding to elements in positions $i$ and $i'$ such that $i$ is linked to $k$ and $i'$ to $k'$.

The variance of $\hat{t}_{y1}$ is given in (2.13), and the $\mathbf{\Delta}$ property states that, for any units $k$ and $k'$ in $U_T$, we have $\mathbf{\Delta}_{kk'} = c_{kk'} \mathbf{1}_k^t \mathbf{1}_{k'}$. By using also the fact that, for all $k \in U_T$, $\tilde{\theta}_k$ satisfies the standardization constraint :

$$
\mathbf{1}_k^t \tilde{\theta}_k = 1 \quad \text{for all} \quad k \in U_T, \tag{2.14}
$$

where $\mathbf{1}_k$ is the $N_{Fk}$-dimensional vector of ones, we have :

$$
\begin{align*}
\mathrm{Var}(\hat{t}_{y1}) &= \sum_{k \in U_T} y_k^2 \tilde{\theta}_k^t \Delta_k \tilde{\theta}_k + \sum_{k \in U_T} \sum_{k' \neq k \in U_T} y_k y_{k'} \tilde{\theta}_k^t \Delta_{kk'} \tilde{\theta}_{k'} \\
&= \sum_{k \in U_T} y_k^2 \tilde{\theta}_k^t \Delta_k \tilde{\theta}_k + \sum_{k \in U_T} \sum_{k' \neq k \in U_T} y_k y_{k'} c_{kk'} \tilde{\theta}_k^t \mathbf{1}_k \mathbf{1}_{k'}^t \tilde{\theta}_{k'} \\
&= \sum_{k \in U_T} y_k^2 \mathrm{Var}(\hat{t}_{\tilde{\theta}_k}) + \sum_{k \in U_T} \sum_{k' \neq k \in U_T} y_k y_{k'} c_{kk'}. \tag{2.15}
\end{align*}
$$

This is exactly (2.6) and finishes the proof of the first part of Proposition 2.2.1.

For the derivation of the optimal link weights, the proof is similar to the proof given
in Deville and Lavallée (2006) (section 6.2) for deriving the optimal weighted links. In our
situation, the links between the frame population $U_F$ and the target population $U_T$ are of
type MtO. Thus, we don't need to use a factorization step as in Deville and Lavallée (2006).

Our aim is to find the link weights $\tilde{\theta}_k^{opt}$, $k \in U_T$ that minimize

$$
\mathrm{Var}(\hat{t}_{y1}) = \sum_{k \in U_T} y_k^2 \tilde{\theta}_k^t \Delta_k \tilde{\theta}_k + \sum_{k \in U_T} \sum_{k' \neq k \in U_T} y_k y_{k'} c_{kk'}
$$

under the standardization constraint (2.14). Since $\Delta_k$, $k \in U_T$ is invertible, this problem is
the minimisation of a positive quadratic form under linear constraint. Thus a solution exists
and is unique. The variance Var($\hat{t}_{y1}$) is minimized for vectors $\tilde{\theta}_k$ verifying the following equation (see Deville and Lavallée, 2006, equation 6.4) :

$$y_k^2 \Delta_k \tilde{\theta}_k = \lambda_k \mathbf{1}_k, \quad k \in U_T, \qquad (2.16)$$

where $\lambda_k, k \in U_T$ the Lagrange multipliers. Let us show that the optimal weights are given by :

$$\tilde{\theta}_k^{opt} = (\mathbf{1}_k^t \Delta_k^{-1} \mathbf{1}_k)^{-1} \Delta_k^{-1} \mathbf{1}_k, \quad k \in U_T,$$

where $\Delta_k = \Delta_{kk}$. Equation (2.16) implies

$$\tilde{\theta}_k = y_k^{-2} \Delta_k^{-1} \lambda_k \mathbf{1}_k, \quad k \in U_T. \qquad (2.17)$$

Multiplying by $\mathbf{1}_k^t$ the equation (2.17) and using again the standardization constraints (2.14), we obtain the following expression for the Lagrange multipliers :

$$\lambda_k = y_k^2 (\mathbf{1}_k^t \Delta_k^{-1} \mathbf{1}_k)^{-1}, \quad k \in U_T. \qquad (2.18)$$

Finally, by plugging the expression of $\lambda_k$ given in (2.18) in the expression of $\tilde{\theta}_k$ from (2.17), we obtain :

$$\begin{aligned} \tilde{\theta}_k^{opt} &= y_k^{-2} \Delta_k^{-1} \mathbf{1}_k \lambda_k \\ &= y_k^{-2} \Delta_k^{-1} \mathbf{1}_k y_k^2 (\mathbf{1}_k^t \Delta_k^{-1} \mathbf{1}_k)^{-1} \\ &= \Delta_k^{-1} \mathbf{1}_k (\mathbf{1}_k^t \Delta_k^{-1} \mathbf{1}_k)^{-1}. \end{aligned}$$

This finishes the proof of 2.2.1. □

## Proof of Proposition 2.2.2

Result (2.15) holds for any standardized set of weights. Thus, the result also holds for the set of optimal weights $\tilde{\theta}_{ik}^{opt}$, and we can write :

$$\mathrm{Var}(\hat{t}_{y1}^{opt}) = \sum_{k \in U_T} y_k^2 \mathrm{Var}(\hat{t}_{\tilde{\theta}_k^{opt}}) + \sum_{k \in U_T} \sum_{k' \neq k \in U_T} y_k y_{k'} c_{kk'}. \qquad (2.19)$$

Using equations (2.15) and (2.19), we get :

$$\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}) - \mathrm{Var}(\hat{t}_{\tilde{\theta}_k^{opt}})).$$