Les différentes transformations des données

add

\bullet Introduction:

Certains outils statistiques nécessitent une transformation des données afin de pouvoir être utilisés ou optimisés. Par exemple, dans le cadre de l’utilisation d’un test de Student ou bien d’une régression linéaire l’hypothèse de normalité est primordiale ou encore dans le cas où les données ne répondent pas aux hypothèses recquises, les équivalents non paramétriques de la théorie des tests nécessitent une transformation des données en vecteur des rangs associés aux données. D’Autres outils encore comme l’analyse en composantes principales, les K-plus proches voisins pondérés ou encore les support vectors machines présentent une certaine sensibilité aux trop grandes dispersions des données à laquelle le statisticien peut pallier.

Les méthodes de transformation des données se déclinent donc en trois familles:

– la standardisation ou l’action de centrer-réduire les données pour diminuer l’échelle de dispersion tout en conservant la forme des distributions conjointes,

– la transformation en vecteur de rangs dont l’objectif est de se concentrer sur l’ordre des valeurs des données et plus sur les valeurs elles-mêmes,

– la normalisation des données qui consistent à rechercher la transformation adaptée au travers de laquelle elles suivront une loi normale.

\bullet Les méthodes classiques:

La standardisation

Il s’agit de centrer-réduire la matrice \mathbf{X} en retranchant pour chaque vecteur X ^j, j \in [1, P] sa moyenne et en divisant par son écart-type. La formule d’usage est alors,

\forall j \in [1, P], X_{Std} ^j = \frac{X ^j - \mu_j}{\sigma_j ^2}

La transformation en rangs

La théorie des tests se divisent en deux familles: les tests paramétriques et les tests non paramétriques. Dans le premier cas X ^j, variable aléatoire continue, doit répondre à des hypothèses d’utilisation sur sa distribution tandis que dans le second cas nous en faisons abstraction. Ainsi, lorsque les hypothèses d’utilisation ne sont pas respectées, les alternatives non paramétriques sont souvent des choix pertinents. Il est alors judicieux de se baser sur une transformation en vecteur de rangs (R ^1, \ldots, R^p) associés aux variables quantitatives continues (X ^1, \ldots, X^p). L’algorithme ci-dessous permet d’expliciter cette transformation pour une variable aléatoire X ^j:

1) ranger les valeurs de X ^j par ordre croissant,

2) chaque R_i ^j correspond à l’indice du rang de X_i ^j dans ce classement,

3) dans le cas où X_i ^j = X_k ^j = \ldots = X_l ^j nous avons:

R_i ^j = R_k ^j = \ldots = R_l ^j = \frac{\sum_{l = i} ^l R_l ^j}{\sharp \lbrace \mbox{Valeurs egales} \rbrace}

La normalisation

Normaliser les données revient à appliquer une fonction f sur X ^j tel que f(X ^j) suit une loi normale. Plusieurs méthodes simples existent comme,

  • la transformation logarithmique: \forall j \in [1, P], X_{Norm} ^j = log (X ^j)
  • l’élévation au carré: \forall j \in [1, P], X_{Norm} ^j = (X ^j) ^2
  • la mise sous racine carré: \forall j \in [1, P], X_{Norm} ^j = \sqrt{X ^j}
  • la transformation arc-sinus: \forall j \in [1, P], X_{Norm} ^j = arcsin \sqrt{X ^j}
  • la mise à l’exponentielle: \forall j \in [1, P], X_{Norm} ^j = e ^{X ^j}
  • la transformation inverse:\forall j \in [1, P], X_{Norm} ^j = \frac{1}{X ^j}

Néanmoins ces méthodes requièrent l’absence de valeurs négatives (à l’exception de la mise à l’exponentielle et de la transformation inverse), l’absence de valeurs nulles (pour la transformation logarithmique et la transformation inverse) ou d’un intervalle spécifique dans lequel X ^j prend ses valeurs (ainsi pour la transformation arc-sinus il faut que X ^j \in [0, 1]).

Dans le cas de la présence de valeurs négatives, un paramètre c de translation peut être appliqué afin de déplacer la distribution sur un intervalle de valeurs positives tout en conservant la forme de la distribution. Pour le cas des valeurs nulles, la méthode à suivre est souvent de « tricher » en supposant qu’il s’agit d’une valeur très faible arrondie et en lui attribuant une valeur décimale (par exemple: 0.000\cdots01) en fonction de l’intervalle de valeurs prises par X ^j de manière à ce qu’elle n’influe pas la distribution générale.

Deux autres outils de normalisation existent et sont plus performants: la transformation Box-Cox et la transformation de Johnson. Nous présentons dans la suite de cette article le premier.

\bullet La normalisation de Box-Cox

C’est en 1964 que George Edward Pelham Box et David Roxbee Cox ont mis à notre disposition ce fabuleux outil qu’est la normalisation (également appelée transformation) de Box-Cox. La transformation de Box-Cox ne permet pas de pallier au problème des valeurs négatives ou nulles pour X ^j. Toutefois, elle offre l’avantage de pouvoir prendre en compte une variable Y qualitative et de normaliser X ^j en fonction des différents groupes de Y. La transformation de Box-Cox est la plus populaire des méthodes de normalisation des données, elle le doit principalement à sa simplicité de mise en oeuvre malgré qu’il puisse exister des situations pour lesquels elle soit inefficace.

La transformation de Box-Cox consiste en réalité en l’application d’un algorithme qui va chercher à maximiser la log-vraisemblance d’une configuration de la variable X ^j en fonction du paramètre de transformation noté \lambda. Il convient de fixer au préalable un intervalle de variation pour \lambda, sachant qu’il n’y a pas de règle encore définie sur ce paramètrage.

L’objectif est donc de maximiser la log-vraisemblance, dépendante de \lambda, et de formule:

LL_{\lambda} = - \frac{n}{2} \cdot log(2 \pi) - n \cdot log(s_{\lambda} ^j) - \sum_{i = 1} ^n [-\frac{1}{2 \cdot s_{\lambda} ^j} \cdot log [ f_{\lambda} (X ^j) - \mu_{\lambda} ^j ] ^2 + (\lambda - 1) \sum_{i = 1} ^n log(X_i ^j)]

\mu_{\lambda} ^j, s_{\lambda} ^j représentent, respectivement, la moyenne et l’écart-type calculés à partir de f_{\lambda} (X ^j) la fonction de transformation de Box-Cox pour \lambda variant dans son intervalle de variation fixé et de formule:

X_{Norm_{\lambda}} ^j = f_{\lambda} (X ^j) = \left\{ \begin{array}{ll} \frac{(X ^j) ^ {\lambda} - 1}{\lambda} & \mbox{si } \lambda \neq 0 \\ log(X ^j) & \mbox{sinon} \\ \end{array} \right.

Notons que le retour au vecteur X ^j se fait via la transformation suivante:

X ^j = (1 + f_{\lambda} (X ^j) \cdot \lambda) ^{\frac{1}{\lambda}}, \forall \lambda \neq 0

Démonstration:

La formule de la log-vraisemblance que nous avons présentée ci-dessus se retrouve assez logiquement. En effet, la densité de probabilité pour X, et donc la vraisemblance en relation avec f_{\lambda}(X), est obtenue en multipliant la densité d’une loi normal aux points de f_{\lambda}(X) avec la jacobienne J de la transformation dont la particularité est:

J(\lambda; X ^j) = f_{\lambda} ' (X) = \frac{\partial (\frac{X ^{\lambda} - 1}{\lambda})}{\partial X} = \frac{\lambda}{\lambda} \cdot X ^{\lambda - 1} = X ^{\lambda - 1}

Par conséquent, si nous partons du produit évoqué ci-dessus avec g(.) fonction de densité d’une loi normale, nous avons:

g(f_{\lambda}(X ^j)) \cdot J(\lambda; X) \Rightarrow L = \prod_{i = 1} ^n g(f_{\lambda} (X_i ^j)) \cdot X_i ^{\lambda - 1}

\Rightarrow LL = log(L) = log \prod_{i = 1} ^n [\frac{1}{s_{\lambda} ^j \sqrt{2 \pi}} e ^{- \frac{1}{2} (\frac{f_{\lambda} (X_i ^j) - \mu_{\lambda} ^j}{s_{\lambda} ^j}) ^2} \cdot X_i ^{\lambda - 1}]

= \sum_{i = 1} ^n log [ \frac{1}{s_{\lambda} ^j \sqrt{2 \pi}} e ^{-\frac{1}{2} (\frac{f_{\lambda} (X_i ^j) - \mu_{\lambda} ^j}{s_{\lambda} ^j}) ^2} X_i ^{\lambda - 1} ]

= \sum_{i = 1} ^n [log(1) - log(s_{\lambda} ^j) - log(\sqrt{2 \pi}) + log( e^ {-\frac{1}{2} (\frac{f_{\lambda} (X_i ^j) - \mu_{\lambda} ^j}{s_{\lambda} ^j}) ^2}) + log(X_i ^{\lambda - 1})]

= \sum_{i = 1} ^n [0 - log(s_{\lambda} ^j) - \frac{1}{2} \cdot log(2 \pi) - \frac{1}{2} (\frac{f_{\lambda} (X_i ^j) - \mu_{\lambda} ^j}{s_{\lambda} ^j}) ^2 + (\lambda - 1) \cdot log(X_i)]

= -n \cdot log(s_{\lambda} ^j) - \frac{n}{2} \cdot log(2 \pi) - \frac{1}{s_{\lambda} ^2} \sum_{i = 1} ^n [ f_{\lambda} (X_i ^j) - \mu_{\lambda} ^j ] ^2 + (\lambda - 1) \cdot \sum_{i = 1} ^n log(X_i)

Autres approches:

Deux alternatives existent, l’une est en fait la version à deux paramètres \lambda_1, \lambda_2 de la transformation de Box-Cox,

X_{Norm} ^j = \left\{ \begin{array}{ll} \frac{(X ^j + \lambda_2) ^ {\lambda_1} - 1}{\lambda_1} & \mbox{si } \lambda_1 \neq 0 \\ log(X ^j + \lambda_2) & \mbox{sinon} \\ \end{array} \right.

La seconde alternative est celle de Johnson, moins populaire et plus complexe, et qui permet notamment de composer lorsque X ^j contient des valeurs négatives ou nulles.

\bullet Exemple:

Soit l’échantillon suivant:

add

Dans un premier temps, nous appliquons les méthodes simples: la transformation logarithmique, l’élévation au carré, la mise sous racine carré, la transformation inverse et la mise à l’exponentielle. Nous obtenons les distributions suivantes:

add2

Nous voyons que les méthodes simples sont assez inefficaces pour normaliser X puisque les différents tests de Shapiro-Wilk effectués rejettent systèmatiquement l’hypothèse de normalité.

Nous proposons désormais d’appliquer l’algorithme de la transformation de Box-Cox. Le but étant de déterminer le \lambda optimisant la normalisation de la variable X. Nous faisons varier notre paramètre dans l’ensemble \lbrace -3, -2, -1, 1, 2, 3 \rbrace.

Nous commençons par la détermination des différents versions normalisées de X pour les différentes valeurs de \lambda ci-dessus via la formule:

f_{\lambda} (X) = \frac{X ^{\lambda} - 1}{\lambda}

, nous obtenons donc la matrice où la première colonne contient les valeurs brutes, les colonnes suivantes les valeurs transformées pour le \lambda qui varie:

\begin{pmatrix} X & \lambda = -3 & \lambda = - 2 & \lambda = - 1 & \lambda = 1 & \lambda = 2 & \lambda = 3 \\ 8.1472 & 0.33271695 & 0.49246726 & 0.87725844 & 7.1472 & 32.68843392 & 179.9285392 \\ 9.0579 & 0.33288480 & 0.49390582 & 0.88959913 & 8.0579 & 40.52277620 & 247.3868031 \\ 1.2699 & 0.17056490 & 0.18995056 & 0.21253642 & 0.2699 & 0.30632300 & 0.3492997 \\ 9.1338 & 0.33289589 & 0.49400669 & 0.89051654 & 8.1338 & 41.21315122 & 253.6663871 \\ 6.3236 & 0.33201512 & 0.48749622 & 0.84186223 & 5.3236 & 19.49395848 & 83.9558639 \\ 0.9754 & -0.02586184 & -0.02553846 & -0.02522042 & -0.0246 & -0.02429742 & -0.0239998 \\ 2.7850 & 0.31790201 & 0.43553565 & 0.64093357 & 1.7850 & 3.37811250 & 6.8670289 \\ 5.4688 & 0.33129534 & 0.48328194 & 0.81714453 & 4.4688 & 14.45388672 & 54.1865438 \\ 9.5751 & 0.33295363 & 0.49454640 & 0.89556245 & 8.5751 & 45.34127001 & 292.2898296 \\ 9.6489 & 0.33296227 & 0.49462950 & 0.89636124 & 8.6489 & 46.05063560 & 299.1082853 \\ 1.5761 & 0.24819455 & 0.29871915 & 0.36552249 & 0.5761 & 0.74204561 & 0.9717254 \\ 9.7059 & 0.33296877 & 0.49469240 & 0.89696988 & 8.7059 & 46.60224740 & 304.4464687 \\ 9.5717 & 0.33295322 & 0.49454252 & 0.89552535 & 8.5717 & 45.30872044 & 291.9782197 \\ 4.8538 & 0.33041837 & 0.47877703 & 0.79397585 & 3.8538 & 11.27968722 & 37.7841639 \\ 8.0028 & 0.33268297 & 0.49219297 & 0.87504373 & 7.0028 & 31.52240392 & 170.5125961 \\ 1.4189 & 0.21664609 & 0.25164871 & 0.29522870 & 0.4189 & 0.50663861 & 0.6188797 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{pmatrix}

Ensuite, calculons pour chaque f_{\lambda}(X), \forall \lambda \in \lbrace -3, -2, -1, 1, 2, 3 \rbrace la moyenne \mu_{\lambda} et l’écart-type s_{\lambda}. Nous pouvons désormais appliquer la formule de la log-vraisemblance en chacun des points f_{\lambda} (X_i), i \in [1, 20] et faisant intervenir les données brutes X_i:

L(f_{\lambda}(X_i)) = -\frac{1}{2 \times s ^2} \times (f_{\lambda} (X_i) - \mu_{\lambda}) ^2 + (\lambda - 1) \times log(X_i)

Ainsi, pour chaque observation et pour chaque valeur du paramètre \lambda, nous obtenons la matrice de résultats suivante:

\begin{pmatrix} LL_{\lambda = -3} & LL_{\lambda = -2} & LL_{\lambda = -1} & LL_{\lambda = 1} & LL_{\lambda = 2} & LL_{\lambda = 3} \\ -8.479135 & -6.409691 & -4.346851 & -0.1369732409 & 2.0182532 & 4.162568432 \\ -8.903788 & -6.732621 & -4.584175 & -0.3195406894 & 1.8673801 & 4.087034199 \\ -1.957951 & -2.142060 & -2.168162 & -1.2184568013 & -0.6638793 & -0.233505522 \\ -8.937220 & -6.758012 & -4.602831 & -0.3381958856 & 1.8445733 & 4.061893685 \\ -7.462282 & -5.632935 & -3.778143 & -0.0004298665 & 1.7958142 & 3.554280961 \\ -6.526091 & -5.158099 & -3.576740 & -1.3617822495 & -0.9517867 & -0.764794082 \\ -4.128943 & -3.075512 & -2.093173 & -0.6070319391 & 0.3298588 & 1.398454525 \\ -6.878034 & -5.184171 & -3.454031 & -0.0415887216 & 1.5293813 & 3.110504307 \\ -9.126233 & -6.901486 & -4.708211 & -0.4571451572 & 1.6766955 & 3.842538392 \\ -9.156986 & -6.924817 & -4.725340 & -0.4787835226 & 1.6424508 & 3.792426287 \\ -1.963800 & -1.776753 & -1.740387 & -1.0778858603 & -0.4166366 & 0.204502335 \\ -9.180577 & -6.942712 & -4.738477 & -0.4958385366 & 1.6147344 & 3.750839250 \\ -9.124810 & -6.900406 & -4.707418 & -0.4561603241 & 1.6782290 & 3.844747729 \\ -6.396877 & -4.813315 & -3.191057 & -0.1127209004 & 1.2960108 & 2.762786138 \\ -8.407442 & -6.355093 & -4.306735 & -0.1150235995 & 2.0234063 & 4.143368587 \\ -1.799064 & -1.825264 & -1.887763 & -1.1489776963 & -0.5385110 & -0.009026059 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \end{pmatrix}

Il ne reste plus qu’à appliquer la formule faisant intervenir la somme de la log-vraisemblance de chaque individu selon les différentes variations du paramètre \lambda considérées:

LL = -\frac{20}{2} \times log(2 \times \pi) - 20 \times log(s_{\lambda}) + \sum_{i = 1} ^{20} L(f_{\lambda} (X_i))

, nous obtenons alors:

\lambda = -3 \Longrightarrow LL = -110.5453,

\lambda = -2 \Longrightarrow LL = -86.86381,

\lambda = -1 \Longrightarrow LL = -67.89594,

\lambda = 1 \Longrightarrow LL = -51.75309,

\lambda = 2 \Longrightarrow LL = -53.52531,

\lambda = 3 \Longrightarrow LL = -58.67262,

Le graphe ci-dessous permet de voir que le paramètre qui maximise la log-vraisemblance est \lambda = 1.

add3

Nous constatons que pour \lambda = 1 nous avons: X_{Norm} = f_{\lambda = 1}(X) = \frac{X ^1 - 1}{1} = X - 1. Nous en déduisons que la normalisation de Box-Cox n’est pas parvenue à trouver une fonction f_{\lambda} permettant à X de suivre une loi normale.

\bullet Application informatique:

Procédure SAS:

– pour la transformation de Box-Cox: http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_transreg_sect015.htm

– pour la transformation de Johnson: http://support.sas.com/documentation/cdl/en/imlsug/63546/HTML/default/viewer.htm#imlsug_ugvartransform_sect002.htm

Package R:

– pour la transformation de Box-Cox: https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/boxcox.html

– pour la transformation de Johnson: https://cran.r-project.org/web/packages/Johnson/index.html

\bullet Bibliographie:

– Transformation de données: normalisation, stabilisation des variances de Daniel Borcard

– Tests de normalité – Techniques empiriques et tests statistiques de  Ricco Rakatomolala

– Estimation of the box-cox transformation parameter and application to hydrologic data de Melanie Wong

– An analysis of transformations de G. E. P. Box et D. R. Cox

– Le document: http://russell.vcharite.univ-mrs.fr/EIE/fchap14.pdf

L’ Imputation des données manquantes

add

\bullet Introduction:

Rare sont les bases de données parfaites et notamment, parfaitement renseignées. Il est régulièrement fort à parier que la collecte des données se fera toujours avec quelques imprévus amenant toujours sa part d’information incomplète. C’est ce que nous appelons les données manquantes ou non-réponse et parfois notée NA (Non Applicable).

La gestion de la non-réponse varie d’un statisticien à l’autre voir d’un domaine d’utilisation à l’autre. Soit elle est conservée dans le jeu de données soit elle est substituée par une estimation. C’est sur ce point précis que la divergence entre utilisateur apparaît. En effet, il peut être philosophiquement encombrant d’estimer une valeur de substitue à la non-réponse et pouvant porter directement atteinte au caractère aléatoire de la distribution.

Les données manquantes peuvent être répertoriées selon trois catégories:

– Missing completely at random (MCAR) : Une donnée est MCAR, c’est à dire manquante de façon complètement aléatoire si la probabilité d’absence est la même pour toutes les observations. Cette probabilité ne dépend que des paramètres extérieurs et indépendants à cette variable. A noter que le cas des données MCAR reste tout de même peu courant.

– Missing at random (MAR) : Il arrive souvent que la non-réponse soit lié à un facteur précis et donc qu’elle ne soit pas entièrement aléatoire. Une donnée est MAR si la probabilité d’absence est liée à une ou plusieurs autres variables observées. Il existe alors des méthodes statistiques appropriées qui permettrons d’éviter de biaiser l’analyse.

– Missing not at random (MNAR) : Une donnée est MNAR si la probabilité d’absence dépend de la variable en question. Les données MNAR induisent une perte de précision (inhérente à tout cas de données manquantes) mais aussi un biais qui nécessite le recours à une analyse de sensibilité.

Les méthodes d’imputations peuvent se diviser en trois principales familles: vulgairement, les méthodes simples basées sur la suppression de la non-réponse ou sa substitution par la moyenne, la médiane, la modalité la plus fréquente ou la création d’une modalité qui lui est propre. Les méthodes déterministes qui consiste à simuler un modèle et à prédire sa valeur à partir de ce dernier et dont le principal défaut est que deux observations aux caractéristiques explicatives identiques donneront la même valeur de substitution. Les méthodes aléatoires qui visent à simuler une valeur de substitution en partant d’une postulat que la distribution suit une structure bien particulière.

\bullet Les méthodes simples:

a. La suppression de l’observation contenant une non-réponse

La méthode est la plus radicale et régulièrement prisée. Elle se résume à supprimer l’observation si elle contient au moins une non-réponse à l’une des caractéristiques considérées dans l’analyse.

Pour la suite des analyses, deux approches sont à prendre en compte:

– soit nous estimons que les résultats doivent être standardisés pour un même effectif total et alors il convient de supprimer une observation si elle a au moins une donnée manquante dans l’une des caractéristiques retenues,

– soit supprimer les observations en fonction des caractéristiques en cours de gestion, comme par exemple pour un modèle de prédiction où nous ne considérerons que les données manquantes des variables du modèle et pas celles rejetées.

b. Transformer la variable

Il s’agit de reconsidérer l’analyse sous un angle différent. La variable se voit attribuer une nouvelle information en codant la non-réponse et en la prenant en compte dans les analyses. Ce type de manipulation est à utiliser uniquement si le taux de non-réponse est important, sinon quoi nous nous retrouvons en présence d’une modalité dont l’effectif est très faible et qui ne peut être prise en compte dans les analyses.

En fonction du format de la variable l’approche s’adapte:

– pour une variable qualitative, une nouvelle modalité qui constituera à représenter la non-réponse,

– pour une variable quantitative, coder la variable en plusieurs modalités dont une qui représente la non-réponse.

c. Imputation par la moyenne, la médiane et la valeur la plus fréquente

Il s’agit d’approches univariées, se basant uniquement sur la variable considérée sans leur interaction avec les autres variables du jeu de données.

Pour des variables quantitatives, l’approche classique est de la remplacer par la moyenne de la distribution si la taille d’échantillon est particulière importante, par la médiane sinon. Pour les variables qualitatives, il s’agit de l’imputer par la modalité la plus prisée.

De telles approches sont particulièrement gênantes étant donné qu’en fonction du taux de non-réponse elles vont modifier considérablement la distribution de la variable et donc les intervalles de confiance associé.

La manipulation peut être améliorée en procédant à une imputation conditionnelle, c’est à dire en prenant en compte une ou plusieurs informations auxiliaires pour adapter la valeur de substitution.

\bullet Les méthodes déterministes

a. Par outil décisionnel

Dans un premier temps, il faut définir un sous-ensemble de variables pour lequel la non-réponse ne coïncide pas avec celle de la variable à imputer et qui soit lié à cette même variable. De même, il ne faut pas que l’observation à imputer ait une caractéristique retenue pour laquelle elle ne soit pas renseignée.

L’idée est simplement de construire un modèle (régression linéaire, logistique, k-plus proches voisins, arbre de décision, etc) et d’en prédire la valeur de substitution à partir de ces caractéristiques renseignées pour l’information manquante à remplacer.

Cette approche présente l’avantage d’être facile à exécuter mais également s’adapte à n’importe quel format de variables.

b. Algorithme LOESS

La régrESSion LOcale consiste à construire un polynôme de degré faible par méthode des moindres carrés pondérés en fonction de la proximité entre les observations jugées semblables à celle sur laquelle nous voulons procéder à l’imputation. La façon dont est conçu l’algorithme implique que cette méthode est à réserver aux variables continues.

Pour l’imputation d’un observation, possédant n_{NR} donnée(s) manquante(s) parmi ses P caractéristiques, à partir de K (paramètre fixé au préalable) voisins, l’algorithme LOESS est,

– Etape N° 1: Chercher les K plus proches voisins.

– Etape N°2: Construire la matrice \mathbf{A} de taille K \times (P - n_{NR}) et tel qu’en ligne nous retrouvons nos K plus proches voisins décrit, en colonne, au travers des P - n_{NR} caractéristiques pour lesquelles notre observation à imputer est renseignée.

– Etape N°3: Construire la matrice \mathbf{B} de taille K \times n_{NR} et qui correspond en réalité à la contraposé de la matrice \mathbf{B}. Elle contient alors en ligne nos K plus proches voisins décrit, en colonne, au travers des n_{NR} caractéristiques pour lesquelles notre observation à imputer n’est pas renseignée.

– Etape N°4: Construire le vecteur P qui correspond aux P - n_{NR} caractéristiques renseignées de l’observation à imputer.

– Etape N°5: Le vecteur u des données manquantes pour cette observation s’écrit alors,

u = \mathbf{B} ^t \cdot (\mathbf{A} ^t) ^{-1} \cdot P,

, où l’inverse de \mathbf{A} s’obtient par pseudo-inversion (car la matrice est assez logiquement rarement carré dans cette situation).

c. L’Algorithme NIPALS

L’Algorithme NIPALS est l’acronyme de Nonlinear Iterative Partial Least Square, faisant référence à la régression PLS et basé sur l’ACP. Grosso-modo, l’idée de cet algorithme est de quadriller l’observation à imputer par les valeurs les plus probables des autres observations. L’approche reste linéaire, puisque basée sur l’ACP, est reste performante uniquement en présence de multicolinéarité entre les variables considérées. De plus, les variables considérées doivent être continues.

En fixant préalablement le nombre de composantes à calculer H, l’algorithme NIPALS se décrit de la manière suivante,

– Etape N°1: Initialiser \mathbf{X} ^0 = \mathbf{X} que nous centrons-réduisons.

– Etape N°2: \forall h \in [1,H]:

————— Initialiser t(h) = (X ^1) ^{h - 1}.

————— Initialiser p(0) = (\frac{1}{\sqrt{P}}, \cdots, \frac{1}{\sqrt{P}}) de taille 1 \times P.

————— Et jusqu’à convergence de p(h),

—————————– Calculer,

p(h) = (\frac{\sum_{i_{renseigne} \in [1,n]} (X_i ^1) ^{h-1} t_i (h)}{\sum_{i_{renseigne} \in [1,n]} t_i (h) ^2}, \cdots, \frac{\sum_{i_{renseigne} \in [1,n]} (X_i ^P) ^{h-1} t_i (h)}{\sum_{i_{renseigne} \in [1,n]} t_i (h) ^2})

, de taille 1 \times P.

—————————— Normaliser p(h) par la formule,

p (h) ^N = (\frac{p_1 (h)}{\sqrt{\sum_{j = 1} ^P p_j (h) ^2}}, \cdots, \frac{p_P (h)}{\sqrt{\sum_{j = 1} ^P p_j (h) ^2}})

—————————— Calculer,

t(h) = (\frac{\sum_{j_{renseigne} \in [1,P]} (X_1 ^j) ^{h-1} p_j(h) ^N}{\sum_{j_{renseigne} \in [1,P]} (p_j (h) ^N) ^2}, \cdots, \frac{\sum_{j_{renseigne} \in [1,P]} (X_n ^j) ^{h-1} p_j (h) ^N}{\sum_{j_{renseigne} \in [1,P]} (p_j (h) ^N) ^2})

, de taille 1 \times n.

—————————— Contrôle de la convergence de p(h) ^N par la formule,

\sum_{j = 1} ^P [p_j(h) ^N - p_j(0) ^N] ^2

—————————— Si il n’y a pas convergence alors nous relançons l’étape en conservant t(h) pour le recalcul de p(h) et l’écrasons par la nouvelle valeur obtenue en fin d’itération. Pour le contrôle de la convergence, p(0) est écrasé par p(h) ^N pour la nouvelle itération. Sinon, nous passons à la phase suivante.

————— Calcul de la matrice résiduelle, \mathbf{X} ^h = \mathbf{X} ^{h - 1} - t (h) \cdot (p (h) ^N) ^t

– Etape N°3: Le calcul des composantes t(h) et des vecteurs de valeurs propres p(h) ^N, \forall h \in [1,H] étant accompli, nous pouvons alors imputer par la valeur,

\sum_{h = 1} ^H t_i(h) p_j(h) ^N

, l’observation i \in [1,n] de la variable j \in [1,P] en n’oubliant pas de « dé-centrer-réduire » les données pour retrouver les valeurs d’origine de \mathbf{X}.

\bullet Les méthode aléatoires

a. Le hotdeck

L’idée du hotdeck est de définir un ensemble d’observations candidates (appelées « donneurs ») et d’imputer la non-réponse en fonction de celle prise par le premier « donneur » trouvé. Il s’agit de considérer un sous-ensemble de variables auxiliaires suffisamment corrélées à la variable d’intérêt pour laquelle l’imputation doit avoir lieu. Nous sélectionnons alors les individus ayant les même caractéristiques auxiliaires et imputons.

Si les variables auxiliaires sont des variables qualitatives ayant un faible nombre de modalités, il est fort probable que plusieurs candidats soient retenus et que leur valeur pour la variable à imputer diverge. Il faut alors choisir un donneur par la liste retenue, ce choix peut se faire selon critère,

– en fonction du donneur le plus proche selon une autre caractéristique (par exemple, régulièrement dans les enquêtes sociales c’est l’adresse d’habitation et sa proximité avec celle de l’observation à imputer qui sera déterminante),

– par tirage aléatoire au sein des donneurs, ce qui implique que la valeur d’imputation varie d’une simulation à l’autre.

Le hotdeck est très prisé car il a l’avantage de ne pas déformer la dispersion de la variable à imputer. De plus il reste très simple à mettre en œuvre, et peut même être étendu à l’imputation de plusieurs valeurs non-renseignées pour une même observation.

Un autre algorithme d’imputation similaire au hotdeck existe: le colddeck. Ce dernier  consiste à chercher des donneurs depuis un autre fichier de données.

b. Inférence bayésienne

Les méthodes basées sur l’inférence bayésienne se reposent sur la connaissance à priori et à postériori de la loi de distribution de X, P(X | \theta), P(\theta | X) avec \theta un échantillon distribué selon les paramètres de la distribution de X.

Globalement, la méthode consiste à,

– fixer la loi de distribution de X,

– estimer les paramètres de la loi fixée selon X,

– pour un grand nombre d’itérations, tirer un échantillon \theta selon la loi de distribution de X et mettre à jour les paramètres avant de renouveller cette étape,

– l’échantillon ainsi construit et suffisamment grand permet alors de comparer la loi théorique avec la loi observée et de fixer la valeur la plus probable par laquelle imputer afin que le modèle observé soit cohérent avec le modèle théorique.

Plusieurs variantes existent et sont basées sur l’inférence bayésienne: l’imputation par méthode de Monte-Carlo par Chaîne de Markov, l’algorithme Esperance-Maximisation (EM), l’algorithme de Gibbs, la méthode de Data Augmentation ou encore l’algorithme EM via approche bootstrap du programme AMELIAII.

\bullet Exemple:

Soit l’échantillon E ci-dessous,

add

Si nous prenons les variables indépendamment, nous obtenons les statistiques descriptives suivantes pour X ^1, X ^2,

\overline{X ^1} = 10.37671, mediane(X ^1) = 10.5575 et sd_{X ^1} = 5.672913

\sharp \lbrace X ^2 = "A" \rbrace = 9 et \sharp \lbrace X ^2 = "B" \rbrace = 9

1. Suppression de toutes les données manquantes

Procédons par la méthode la plus simple, la suppression des données manquantes relatives à notre jeu de données. Les observations N°3 et 15 ont une donnée manquante pour la variable X ^1 et les observations N°4 et 17 pour la variable X ^2.

Après suppression de ces quatre observations nous obtenons les statistiques descriptives suivantes pour les variables X ^1 et X ^2 et les 16 observations retenues,

\overline{X ^1} = 10.99134, mediane(X ^1) = 11.5403 et sd_{X ^1} = 5.70679

\sharp \lbrace X ^2 = "A" \rbrace = 8 et \sharp \lbrace X ^2 = "B" \rbrace = 8

2. Imputation par la moyenne, médiane, valeur la plus fréquente

Dans un premier temps procédons à une imputation par la moyenne et ensuite la médiane pour X ^1. Nous avons déjà ces valeurs dans la présentation de notre échantillon et pouvons rester sur un échantillon de 20 observations par conséquent.

Les nouvelles statistiques descriptives après imputation par la moyenne,

\overline{X ^1} = 10.37671, mediane(X ^1) = 10.37671 et sd_{X ^1} = 5.366039

, et après imputation par la médiane,

\overline{X ^1} = 10.39479, mediane(X ^1) = 10.39479 et sd_{X ^1} = 5.366327

Pour X^2, nous remarquons que les modalités « A » et « B » ont la même fréquence d’apparition. Aléatoirement nous optons pour une imputation par la modalité « B » et obtenons alors les statistiques descriptives suivantes,

\sharp \lbrace X ^2 = "A" \rbrace = 8 et \sharp \lbrace X ^2 = "B" \rbrace = 10

3. Par algorithme décisionnelle

Pour l’imputation des valeurs manquantes de X ^2 nous allons privilégier la régression logistique de X ^2 sur X ^1. Le modèle construit sur les 16 observations est le suivant,

P(X ^2 = "B") = \frac{e ^{-0.45690 + 0.04156 \times X ^1}}{1 + e ^{-0.45690 + 0.04156 \times X ^1}}

Imputons désormais les observations manquantes N°4 et 17. Nous avons alors,

P(X_4 ^2 = "B") = \frac{e ^{-0.45690 + 0.04156 \times X_4 ^1}}{1 + e ^{-0.45690 + 0.04156 \times X_4 ^1}}

= \frac{e ^{-0.45690 + 0.04156 \times 4.0135}}{1 + e ^{-0.45690 + 0.04156 \times 4.0135 ^1}}

= \frac{0.7481895}{1.748195}

= 0.4279783 < 0.5

\Rightarrow X_4 ^2 = "A"

P(X_{17} ^2 = "B") = \frac{e ^{-0.45690 + 0.04156 \times X_{17} ^1}}{1 + e ^{-0.45690 + 0.04156 \times X_{17} ^1}}

= \frac{e ^{-0.45690 + 0.04156 \times 6.9059}}{1 + e ^{-0.45690 + 0.04156 \times 6.9059}}

= \frac{0.843757}{1.843757}

= 0.4576292 < 0.5

\Rightarrow X_{17} ^2 = "A"

Nous obtenons alors les nouvelles statistiques descriptives pour la variable X ^2 et sur nos 18 observations,

\sharp \lbrace X ^2 = "A" \rbrace = 10 et \sharp \lbrace X ^2 = "B" \rbrace = 8

Pour l’imputation des valeurs manquantes de X ^1 nous allons privilégier la régression linéaire pour variable explicative qualitative de X ^2 sur X ^1. Le modèle construit sur les 16 observations est le suivant,

Y = 10.363 + 1.257 \times 1_{X ^2 = B}

Imputons désormais les observations manquantes N°3 et 15. Nous avons alors,

X_3 ^1 = 10.363 + 1.257 \times 1_{X ^2 = B} = 10.363 + 1.257 \times 1 = 11.62

X_{15} ^1 = 10.363 + 1.257 \times 1_{X ^2 = B} = 10.363 + 1.257 \times 0 = 10.363

Nous obtenons alors les nouvelles statistiques descriptives pour la variable X ^1 et sur nos 18 observations,

\overline{X ^1} = 10.43819, mediane(X ^1) = 10.7656 et sd_{X ^1} = 5.373245

4. Algorithme LOESS

Nous allons maintenant imputer les données manquantes pour les observations 3 et 15 pour la variable X ^1 en nous aidant de la variable X ^3. Fixons K = 2 voisins sur lesquels nous appuyer.

Les deux voisins les plus proches de l’observation 3 sont les observations 14 et 5. Nous avons alors,

A = (1.935572, 2.115756)

, qui correspond aux valeurs de la variable X ^3 pour les observations 5, 14 et seule variable continue pour laquelle l’observation à imputer est renseignée.

Et,

B = (5.0515, 13.9288)

, qui correspond aux valeurs de la variable X ^1 pour les observations 5, 14 et seule variable pour laquelle l’observation à imputer est non renseignée.

Il nous reste plus qu’à déterminer le dernier terme, et qui correspond à la valeur de l’observation à imputer pour X ^3, car seule variable continue pour laquelle elle est renseignée,

P = 0.2839578

Maintenant que nous avons nos trois éléments, il ne reste plus qu’à calculer la valeur d’imputation au travers de la formule,

u = (5.0515,13.9288) \times (\begin{pmatrix} 1.935572 \\ 2.115756 \\ \end{pmatrix}) ^{-1} \times 0.2839578

= (5.0515,13.9288) \times (\begin{pmatrix} 0.2353891 \\ 0.2573017 \\ \end{pmatrix}) \times 0.2839578

= 4.772972 \times 0.2839578

= 1.355323

Et donc,

X_3 ^1 = u = 1.355323

En procédant à un calcul similaire pour X_{15} ^2 et dont les deux plus proches voisins sont les observations 5, 14 (oui ce sont les mêmes, simple coïncidence), nous obtenons,

X_{15} ^1 = 0.7991907

Les nouvelles statistiques descriptives sur nos 20 observations sont alors \overline{X ^1} = 9.448591, mediane(X ^1) = 9.4661 et sd_{X ^1} = 6.081304.

5. L’algorithme NIPALS

Nous allons imputer les valeurs manquantes de X ^1 à l’aide de l’algorithme NIPALS. Nous travaillons donc sur le couple \mathbf{X} = (X ^1, X ^3) centré-réduit de variables continues et fixons le nombre de composantes à calculer à deux.

Commençons l’algorithme avec l’étape h = 1. Nous initialisons les différents éléments,

\mathbf{X} ^0 = \mathbf{X}, t(1) = (X ^1) ^0, p(0) = (\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}) = (0.7071068,0.7071068)

Nous avons alors,

p(1) = (\frac{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} (X_i ^1) ^0 \times t_i (1)}{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} t_i (1) \times t_i (1)}, \frac{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} (X_i ^2) ^0 \times t_i (1)}{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} t_i (1) \times t_i (1)})

= (\frac{18.83657}{18.83657},\frac{-2.031762}{18.83657})

= (1,-0.1078627)

, que nous normalisons et obtenons ainsi,

p(1) ^N = (0.9942331,-0.1072406)

Ainsi que,

t(1) = (\frac{\sum_{i \in [1,20]} (X_i) ^0 \times p_i (1) ^N}{\sum_{i \in [1,2]} p_i (1) ^N \times p_i (1) ^N}, \cdots, \frac{\sum_{i \in [1,20]} (X_i ^2) ^0 \times p_i (1) ^N}{\sum_{i \in [1,2]} p_i (1) ^N \times p_i (1) ^N})

= (\frac{-1.852308}{1}, \cdots, \frac{1.645699}{1})

= (-1.852308, \cdots, 1.645699)

Nous vérifions si le vecteur p(1) ^N a déjà convergé,

\sum_{j \in [1,2]} [p_j(1) ^N - p_j(0)] ^2 = 0.2871263 ^2 + (-0.8143474) ^2 = 0.7456032 \neq 0

Nous relançons une seconde étape de calcul pour la même itération h = 1 avec p(0) = p(1)_N. Nous obtenons alors,

p(1) ^N = (0.9766683, -0.2147535)

, et,

t(1) = (-1.90459232, \cdots, 1.65469395)

Nous vérifions à nouveau la convergence,

\sum_{j \in [1,2]} [p_j(1) ^N) - p_j(0)] ^2 = (-0.01756481) ^2 + (-0.10751285) ^2 = 0.0118675 \neq 0

Nous n’avons toujours pas convergé, nous devons donc relancer à nouveau une nouvelle itération de calcul pour l’étape h = 1 avec p(0) = p(1) ^N et t(0) = t(1). Au bout de 42 itérations le vecteur p(1) ^N converge enfin et nous avons alors,

p(1) ^N = (0.5115357, -0.8592620)

, et,

t(1) = (-1.57778207, \cdots, 1.12650181)

Nous pouvons désormais boucler l’étape h = 1 par le calcul de la nouvelle matrice résiduelle,

\mathbf{X} ^1 = \mathbf{X} ^0 - t(1) \times p(1) ^N

= \begin{pmatrix} -1.77924460 & 0.77698647 \\ -1.55313432 & -1.19148868 \\ NA & 0.02313673 \\ -1.19098791 & -0.46961838 \\ \cdots & \cdots \\ -0.64503133 & 0.75615571 \\ 1.47792131 & -1.0698482 \\ 1.63309726 & -1.48294518 \\ 1.61771368 & -0.34795388 \end{pmatrix} - \begin{pmatrix} -0.80709182 & 1.35572824 \\ 0.11730402 & -0.1704372 \\ -0.01377375 & 0.02313673 \\ -0.10522675 & 0.7675669 \\ \cdots & \cdots \\ -0.50114763 & 0.84181252 \\ 0.85696750 & -1.43950790 \\ 1.07914897 & -1.81272155 \\ 0.57624587 & -0.96796025 \end{pmatrix}

= \begin{pmatrix} -0.9721579 & 0.5787418 \\ -1.67043834 & -0.9944450 \\ NA & -0.000000000000000003469447 \\ -1.08576115 & -0.6463751 \\ \cdots & \cdots \\ -0.14388370 & -0.08565681 \\ 0.62095381 & -0.3696661 \\ 0.55394829 & 0.3297764 \\ 1.04146781 & -0.6200064 \end{pmatrix}

Procédons maintenant à la second étape h = 2. En reprenant les éléments calculées pendant la première étape avec cette fois-ci t(2) = (X ^1) ^1. Nous avons,

p(2) = (\frac{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} (X_i ^1) ^1 \times t_i (1)}{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} t_i (1) \times t_i (1)}, \frac{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} (X_i ^2) ^1 \times t_i (1)}{\sum_{i \in [1,20] - \lbrace 3, 15 \rbrace} t_i (1) \times t_i (1)}) = (1,0.5953198)

, que nous normalisons et obtenons ainsi,

p(2) ^N = (0.8592620,0.5115357)

Ainsi que,

t(2) = (\frac{\sum_{i \in [1,20]} (X_i) ^1 \times p_i (2) ^N}{\sum_{i \in [1,2]} p_i (2) ^N \times p_i (2) ^N}, \cdots, \frac{\sum_{i \in [1,20]} (X_i ^2) ^0 \times p_i (2) ^N}{\sum_{i \in [1,2]} p_i (2) ^N \times p_i (2) ^N}) = (-1.131381, \cdots, 1.212049)

Nous vérifions la convergence de p(2) ^N,

\sum_{j \in [1,2]} [p_j (2) ^N - p_j(0)] ^2 = 0.1521553 ^2 + (-0.1955711) ^2 = 0.06139928 \neq 0

Nous devons relancer une nouvelle itération de calcul afin d’obtenir la convergence recherchée, nous fixons alors p(0) = p(2)_N et obtenons,

p(2) ^N = (0.8592620,0.5115357)

Et,

t(2) = (-1.131381, \cdots, 1.212049)

Nous vérifions la convergence de p(2)_N pour cette seconde itération,

\sum_{j \in [1,2]} (p_j (2) ^N - p_j(0)) ^2 = 0^2 + 0 ^2 = 0

, qui a donc enfin convergé.

Nous avons donc nos deux composantes p(1) ^N, p(2) ^N et nos deux vecteurs de valeurs propres t(1), t(2). Il nous reste plus qu’à estimer les valeurs à imputer pour les observations 3, 15 de X ^1,

(X_3 ^1)' = \sum_{h = 1} ^2 t_3 (h) \times p_1 (h) ^N

= -0.02692628 \times 0.5115357 + (-0.000000000000000006782414) \times 0.8592620

= -0.1377375

\Rightarrow X_3 ^1 = (X_3 ^1)' \times 5.2975872 + 10.3231839 = 9.593507 une fois « dé-centré-réduit ».

X_{15} ^1 = \sum_{h = 1} ^2 t_{15} (h) \times p_1 (h) ^N

= 0.12300570 \times 0.5115357 + 0 \times 0.8592620

= 0.06292181

\Rightarrow X_{15} ^1 = (X_{15} ^1)' \times 5.2975872 + 10.3231839 = 10.65652 une fois « dé-centré-réduit ».

Nous obtenons alors les nouvelles statistiques descriptives suivantes pour X ^1 et pour ses 20 observations: \overline{X ^1} = 10.30337, mediane(X ^1) = 10.30166, sd(X ^1) = 5.278193.

6. Le hotdeck

Maintenant procédons à une imputation par hotdeck. Commençons par la variable X ^2 en cherchant un donneur à partir de Y en ordonnant l’échantillon en fonction de X ^1.

– Pour l’observation N° 4, le donneur le plus proche et de même valeur pour Y est l’observation N° 2 de valeur X_2 ^1 = "B". Par conséquent X_4 ^2 = "B".

– Pour l’observation N° 17, qui passe alors en 6ème position, le donneur le plus proche et de même valeur pour Y est l’observation N°5 de valeur X_5 ^2 = "A". Par conséquent X_{17} ^2 = "A".

Nous obtenons alors les nouvelles statistiques descriptives pour la variable X ^2 et sur nos 20 observations,

\sharp \lbrace X ^2 = "A" \rbrace = 9 et \sharp \lbrace X ^2 = "B" \rbrace = 9

Procédons à la même manipulation pour la variable X ^1 mais en conservant l’ordre initial.

– Pour l’observation N° 3, le donneur le plus proche et de même valeur pour Y est l’observation N°2 de valeur X_2 ^1 = 2.0949. Par conséquent X_3 ^1 = 2.0949

Pour l’observation N° 15, le donneur le plus proche et de même valeur pour Y est l’observation N°14 de valeur X_{14} ^1 = 13.9288. Par conséquent X_{15} ^1 = 13.9288

Nous obtenons alors les nouvelles statistiques descriptives pour la variable X ^1 et sur nos 20 observations,

\overline{X ^1} = 10.14022, mediane(X ^1) = 10.5575 et sd_{X ^1} = 5.745388

\bullet Application informatique:

Procédure SAS:

– (pour les méthodes d’imputation) https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mi_sect004.htm

– (pour la combinaison des résultats) http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mianalyze_sect004.htm

Fonction R: http://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-scenar-app-idm.pdf

\bullet Bibliographie

– Comprendre et utiliser les statistiques dans les sciences de la vie de Bruno Falissard

– Data mining et statistique décisionnelle. L’intelligence des données de Stéphane Tufféry

– Les techniques de sondage de Pierre Ardilly

– Probabilité, analyse des données et statistique de Gilbert Saporta

– Le document: http://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-m-app-idm.pdf

– Les données manquantes en bio-statistique de N. Meyer

– Données manquantes ou censurées : principes de base de Bernard Delyon

– Processus stochastique et applications de Nicolas Bouleau

– Données manquantes et prévisions: méthodes à imputation variable de Antonio Anselmi, Paola M. Chiodini et Flavio Verrecchia

– Analyse statistique des données longitudinales de Jean-Jacques Droesbeke et Gilbert Saporta

– The Calculation of Posterior Distributions by Data Augmentation de Martin A. Tanner et Wing Hung Wong

Le redressement des pondérations d’un échantillon

add\bullet Introduction:

Afin de saisir le principe général du redressement des pondérations il faut se placer dans la théorie des sondages. L’idée d’un plan de sondage est de déterminer la stratégie optimale visant à sélectionner une sous-population E de taille n selon certains critères et qui sera représentative de la la population d’étude de taille N. La raison d’une telle stratégie étant simple et évidente: diminuer le coût d’enquête et plus particulièrement celui du déploiement humain. De plus, il a été prouvé qu’une enquête par sondage pouvait s’avérer plus performante qu’un recensement (qui vise à enquêter toute la population d’étude) et pouvant apporter un taux de non-réponse important et sur lequel le sondeur perd indirectement tout contrôle.

Ainsi, une fois l’enquête terminée, une unité statistique de la base de sondage  devient la voix d’un certain nombre d’unités statistiques non sondées et de profils similaires et se voit attribuer un poids de sondage. Le redressement des pondérations vise donc à réajuster ces poids de sondage initiaux et liés au plan de sondage et, selon des critères définis au préalable, assurer la reproductibilité de la structure de la population d’étude sur la sous-population tirée.

Il convient de toujours parler de « redressement des pondérations » et non simplement de « redressement ». En effet, dans le domaine de la statistique et de la probabilité, le redressement peut également évoquer celui de la non-réponse et qui représente une gamme d’outils à part entière.

Le grand principe des méthodes de redressement des pondérations est assez logique. Il s’agit de partir des totaux (appelés marges) observés sur l’échantillon tiré et de trouver les pondérations visant à assurer l’égalité avec les marges de la population d’étude. Une telle approche soulève forcément un semblant de paradoxe: à quoi sert-il de faire un sondage si j’ai déjà l’information que je veux mesurer? La réponse permet de justifier l’aspect indirecte et fondamental du redressement des pondérations. En fait, il s’agit de déterminer une série de variables de calage (appelées également variables auxiliaires dans le jargon des plans de sondage) sur lesquelles nous pouvons avoir accès à l’information facilement pour la population d’étude et qui est également corrélées à la thématique de l’étude. Par exemple, dans plusieurs études sociales, le sondeur utilise souvent comme variables de calage: le sexe, le niveau de diplôme et les tranches d’âge depuis les résultats du recensement de la population fournis par l’INSEE. Ces variables sont souvent considérées comme les plus neutres et les plus liées aux différentes thématiques sociales.

Les différentes principales méthodes de redressement des pondérations sont:

– la post-stratification simple: lorsque nous avons une ou plusieurs variable(s) auxiliaire(s) qualitative(s),

– l’estimation par le quotient: lorsque nous avons une variable auxiliaire quantitative,

– l’estimation par la régression: lorsque nous avons plusieurs variables auxiliaires quantitatives;

– le calage sur marge, soit la généralisation des trois méthodes ci-dessus. 

\bullet Les méthodes classiques:

Nous commencerons par évoquer les deux principales méthodes connues: la post-stratification et l’estimation par la régression (dont la méthode d’estimation par quotient est un cas particulier).

1) La post-stratification

De toutes les méthodes utilisées pour redresser les pondérations, celle par post-stratification est la plus célèbre. L’algorithme (qui porte le nom de Raking Ratio), consiste simplement à déterminer les pondérations qui permettent de retrouver les marges attendues à partir de celles observées sur l’échantillon. Il reste assez trivial à comprendre et à mettre en oeuvre. Pour cela, il faut considérer le tableau croisé dynamique issu des variables qualitatives observées sur l’échantillon et la table des marges attendues: M.

L’algorithme consiste en un ajustement itératif des différentes marges, les unes après les autres, associées aux différentes modalités k_p des différentes p \in [1, P] variables auxiliaires utilisées. Ainsi, nous avons à l’étape t_p le réajustement suivant pour la marge de la modalité k_p de la variable X ^p sachant que nous avons fixés les p-1 autres variables:

\forall k_p \in [1, K_p], n_{X ^1 = k_1, \cdots, X ^p = k_p, \cdots, X ^P = k_P} ^t = n_{X ^1 = k_1, \cdots, X ^p = k_p, \cdots, X ^P = k_P} ^{t _ 1} \times \frac{M_{X ^p = k_p}}{n_{X ^p = k_p}}

Une fois l’ajustement effectué pour t_p, si les marges recalculées ne sont pas égales à celles attendues (M) nous poursuivons pour t_{p+1}. Si p + 1 > P alors nous reprenons le processus pour la première variable a temps t + 1 et noté (t + 1)_1L’algorithme tourne jusqu’à ce que les marges attendues soit retrouvées.

Les pondérations W redressées sont alors affectées en fonction du croisement de modalités, auxquels appartient l’observation que nous voulons pondérer, et des effectifs n ^t en l’état au moment de l’arrêt de l’algorithme. Nous avons donc:

W_{i; X ^1 = k_1, \cdots, X ^P = k_P} = \frac{1}{n} \cdot \frac{n_{X ^1 = k_1, \cdots, X ^P = k_P} ^t}{n_{X ^1 = k_1, \cdots, X ^P = k_P}}

2) L’Estimation par la régression

L’estimation par quotient (également appelé estimateur par ratio) est un cas particulier de l’estimation par la régression qui est adaptée au cas des variables quantitatives.

L’idée est de partir du principe que pour \mathbf{X} matrice de variables auxiliaires observée, \mathbf{Y} variables d’intérêt attendues et \mathbf{U} résidus  dont la somme donne 0 (et qui par convention, sont des valeurs proches de 0), nous avons l’équation suivante qui relie \mathbf{X} à \mathbf{Y}:

\forall i \in E; \mathbf{Y}_i = a + b \cdot \mathbf{X}_i + \mathbf{U}_i

Il existe une infinité de couples (a, b) solutions du système ci-dessus. Par exemple, si nous nous retrouvons dans le cas de l’estimateur par quotient pour une variable auxiliaire, alors le couple solution peut être (0, \frac{\overline{Y}}{\overline{X}}). Une méthode plus générale, celle des moindres carrés, donne la formule type suivante pour déterminer un couple solution:

\hat{\hat{b}} = \frac{\sum_{i = 1} ^n (\mathbf{X}_i - \overline{\mathbf{x}}) \cdot (\mathbf{Y}_i - \overline{\mathbf{y}})}{\sum_{i = 1} ^n (\mathbf{X}_i - \overline{\mathbf{x}}) ^2}

\hat{\hat{a}} = \overline{\mathbf{y}} - \hat{\hat{b}} \cdot \overline{\mathbf{x}}

Par conséquent, les poids de redressement W sont de la forme,

W_i = \frac{1}{n} + (\overline{\mathbf{X}} - \overline{\mathbf{x}}) \cdot \frac{\mathbf{X}_i - \overline{\mathbf{x}}}{\sum_{i \in E} (\mathbf{X}_i - \overline{\mathbf{x}}) ^2}

, qui peuvent être négatifs étant donné la forme de W_i. Cette situation reste tout de même rare et sans réelle logique si nous nous limitons au principe de la pondération. Néanmoins, au sens statistique, la conservation de pondérations négatives trouve toute son importance dans le gain de variance que peut amener la méthode d’estimation par la régression. 

\bullet Le CALage sur MARges:

Le calage sur marge est une généralisation de l’ensemble des méthodes de redressement des pondérations. Nous lui dédions sa propre partie.

Le problème:

Utiliser le calage sur marge revient à trouver une solution au problème suivant:

Soit notre échantillon E de taille n obtenu suite à une méthode de sondage quelconque (ce qui constitue une hypothèse d’utilisation essentiel tout de même) depuis notre population d’étude de taille N. Nous disposons alors des pondérations initiales p_i, \forall i \in E et du tableau des marges (ou pour être plus concret: celui des totaux)  attendues M et correspondant aux variables auxiliaires choisies et sur lesquelles caler notre échantillon. Nous cherchons donc à redresser nos pondérations initiales p = (p_1, \cdots, p_n) vers de nouvelles pondérations W = (W_1, \cdots, W_n) tel que la structure de la population d’étude soit retrouvée à partir de la pondération par W de E.

La contrainte à respecter est alors que les  | W_i - p_i |, \forall i soient les plus proches possibles de 0 étant donné que les p_i ont été déterminées selon un plan de sondage garantissant des estimateurs sans biais.

La procédure:

Nous partons donc du système d’équations, dépendant des poids initiaux p, suivant:

\begin{cases} \sum_{i \in E} p_i \cdot X_i ^1 \neq \sum_{i = 1} ^N X_i ^1 \\ \vdots \\ \sum_{i \in E} p_i \cdot X_i ^P \neq \sum_{i = 1} ^N X_i ^P \end{cases}

, et cherchons à déterminer le vecteur des W tel que,

\begin{cases} \sum_{i \in E} W_i \cdot X_i ^1 = \sum_{i = 1} ^N X_i ^1 \\ \vdots \\ \sum_{i \in E} W_i \cdot X_i ^P = \sum_{i = 1} ^N X_i ^P \end{cases}

Pour cela, nous commençons par introduire une distance D(.,.) reliant les pondérations initiales p_i aux pondérations finales W_i, ce qui nous mène au problème d’optimisation suivant:

– minimiser \sum_{i \in E} D(W_i, p_i),

– sous contrainte que \forall p \in [1,P], \sum_{i \in E} W_i \cdot X_i ^p = \sum_{i = 1} ^N X_i ^p.

Ce problème se résout en deux temps:

1) Fixer une fonction F(.), reliée à D, et ayant quelques bonnes propriétés mathématiques. Nous résolvons alors l’équation suivante en déterminant le vecteur des multiplicateurs de Lagrange \lambda = (\lambda_1, \cdots, \lambda_P) et solution du système:

\begin{cases} \sum_{i \in E} p_i \cdot F(X_i ^1 \cdot \lambda_1) \cdot X_i ^1 = \sum_{i = 1} ^N X_i ^1 \\ \vdots \\ \sum_{i \in E} p_i \cdot F(X_i ^P \cdot \lambda_P) \cdot X_i ^P = \sum_{i = 1} ^N X_i ^P \end{cases}

2) Une fois le vecteur de nombres réels \lambda déterminés, nous pouvons calculer:

\forall i \in E: W_i = p_i \cdot F(\mathbf{X}_i \cdot \lambda)

Notons qu’il n’existe pas vraiment de critère pour le choix de la distance D et de la fonction lien F puisque les méthodes restent équivalentes. Généralement, nous choisissons la méthode pour laquelle nous obtenons:

– le plus de poids finaux W proches des poids initiaux p

– le moins de poids finaux W < 0,

– une dispersion des poids finaux var(W) \rightarrow 0.

Les différentes fonctions de lien F

F(.) est donc assimilable au coefficient multiplicateur des poids d’origine p permettant le redressement en W. Les différentes formes qui relient F à D sont:

– La méthode Linéaire: D(W_i, p_i) = \frac{(W_i - p_i) ^2}{p_i} et F(X) = 1 + X.

Nous retrouvons là le cas particulier de l’estimation par régression. La méthode linéaire présente l’avantage d’être très simple pour la résolution du système, contrairement aux autres méthodes qui opèrent par approximation. A noter qu’il s’agit de la seule méthode pouvant générer des pondérations redressées négatives.

– La méthode du Raking Ratio: D(W_i, p_i) = W_i \cdot Log (\frac{W_i}{p_i}) + p_i - W_i et F(X) = e ^X.

Nous retrouvons là le cas particulier de la post-stratification par algorithme du Raking Ratio. La méthode génère des pondérations exclusivement positives mais ne permet pas de contrôler plus optimalement la dispersion des W que la méthode linéaire.

– La méthode linéaire tronquée pour le couple de bornes (B_{inf}, B_{sup}) \in R ^+ fixé au préalable: D(W_i, p_i) est égale à la méthode linéaire si B_{inf} \leq \frac{W_i}{p_i} \leq B_{sup}, sinon +\infty.

Et F(X) est égale à la méthode linéaire si B_{inf} - 1 \leq X \leq B_{sup} - 1, à B_{inf} si X < B_{inf} - 1, B_{sup} si X > B_{sup} - 1.

– La méthode logistique pour le couple de bornes (B_{inf}, B_{sup}) \in R ^+ fixé au préalable: D(W_i,p_i) = (W_i - B_{inf} \cdot p_i) \cdot log(\frac{\frac{W_i}{p_i} - B_{inf}}{1 - B_{inf}}) + (B_{sup} \cdot p_i - W_i) \cdot log(\frac{B_{sup} - \frac{W_i}{p_i}}{B_{sup} - 1}) si B_{inf} \leq \frac{W_i}{p_i} \leq B_{sup}, sinon +\infty.

Et F(X) = \frac{B_{inf} \cdot (B_{sup} - 1) + B_{sup} \cdot (1 - B_{inf}) e ^{c \cdot X}}{B_{sup} - 1 + (1 - B_{inf}) \cdot e ^{c \cdot X}} avec c = \frac{B_{sup} - B_{inf}}{(1 - B_{inf}) \cdot (B_{sup} - 1)}.

A noter que la définition de la méthode logistique en fait également une méthode tronquée. Ces dernières ont pour avantage de pouvoir fixer les bornes (B_{inf}, B_{sup}) qui limitent la dispersion des pondérations W et empêchent ainsi des pondérations trop petites ou trop grandes.

\bullet Présentation et convergence pour F et D:

Nous proposons ici de rentrer plus en détail dans les quatre différentes approches du calage sur marges. Chaque approche a ses qualités et ses défauts, ce qui est la raison principale de l’absence de réel critère de sélection.

Soit,

\mathbf{X} la matrice composée de vecteurs de variables continues et/ou qualitatives. En présence de variables qualitatives il faudra considérer dans \mathbf{X} le tableau disjonctif complet issu de ces variables,

p = (p_1, \cdots, p_n) le vecteur des poids initiaux associés aux différentes unités statistiques de l’échantillon E, soit ceux calculés suite au plan de sondage,

M le vecteur contenant les marges de la population d’étude. Dans le cadre de variables qualitatives de calage, il s’agira des effectifs pour les différentes modalités par variable qualitative. Dans le cadre de variable quantitative, il s’agira de leur somme respective.

W = (W_1, \cdots, W_n) le vecteur des poids redressés finaux.

Comme nous l’avons vu, l’objectif du calage sur marges est de déterminer les multiplicateurs de Lagrange \lambda et solution du système:

\begin{cases} \sum_{i \in E} p_i \cdot F(X_i ^1 \cdot \lambda_1) \cdot X_i ^1 = \sum_{i = 1} ^N X_i ^1 \\ \vdots \\ \sum_{i \in E} p_i \cdot F(X_i ^P \cdot \lambda_P) \cdot X_i ^P = \sum_{i = 1} ^N X_i ^P \end{cases}

, chaque méthode dispose de son propre algorithme afin d’y parvenir.

1) Pour la méthode linéaire:

Soit,

\lambda = [(\mathbf{X} \cdot p) ^t \cdot \mathbf{X}] ^{-1} \cdot [\mathbf{M} - p ^t \cdot \mathbf{X}]

, par conséquent le vecteur des poids finaux est:

W = 1 + (\mathbf{X} \cdot \lambda)

Comme nous l’avions précisé, la méthode linéaire est simple. Les poids peuvent en effet être négatifs car \lambda dépend du terme \mathbf{M} - p ^t \cdot \mathbf{X} qui peut être négatif si la marge attendue est inférieur à celle observée. Ce résultat est assez logique puisque l’idée du redressement des pondérations est de gonfler ces dernières pour faire tendre les totaux vers une valeur plus grande (N >>>> n). Dans le cas où la valeur attendue est plus petite que la valeur observée il reste assez logique que le modèle va intégrer des termes négatifs afin d’équilibrer l’équation linéaire.

Concernant les éventuels problèmes de convergence, l’algorithme se base sur l’estimation des paramètres d’un modèle linéaire, elle est donc assurée quelque soit la situation.

2) La méthode du Raking Ratio:

La méthode du Raking Ratio est un algorithme itératif de résolution de \lambda.

La phase d’initialisation consiste à mettre \lambda, W = 0, fixer le nombre d’itérations maximale maxI et définir la tolérance \epsilon qui servira à déterminer si les marges recalculées sont suffisamment proches des marges attendues pour arrêter l’algorithme et décider s’il y a eu convergence.

Donc, à l’itération t \leq maxI, nous avons tout d’abord la procédure de vérification visant à comparer les marges recalculées avec les pondérations W ^{t - 1} et les marges attendues de M. Si elles ne sont pas toutes égales à \epsilon prés nous déployons la procédure ci-dessous pour l’étape en cours:

\begin{cases} \mbox{a - Calculer: } \lambda ^t = \lambda ^{t - 1} - [(\mathbf{X} \cdot W ^{t - 1})' \cdot \mathbf{X}] ^{-1} \cdot (\mathbf{X}' \cdot W ^{t - 1} - M) \\ \mbox{b - Calculer: } W ^t = p \cdot e ^{\mathbf{X} \cdot \lambda ^t} \end{cases}

L’algorithme Raking Ratio génère des poids exclusivement positif étant donné le passage à l’exponentielle pour la détermination de W.

L’algorithme diverge si \nexists t tel que | max(\frac{\mathbf{X} ' \cdot W ^t - \mathbf{M}}{\mathbf{M}}) | < \epsilon. Il suffit donc qu’une seule marge recalculée soit trop éloignée de celle attendue. Or, il s’agit là probablement du principal défaut de cet algorithme, l’application de la fonction exponentielle aura tendance à faire exploser un peu trop rapidement le vecteur W l’empêchant de converger.

3) La méthode linéaire tronquée:

La méthode linéaire tronquée est un algorithme itératif de résolution de \lambda.

La phase d’initialisation consiste à fixer les bornes B_{inf}, B_{sup} de l’intervalle de variation de W, puis à initialiser \lambda, W avec la procédure de la méthode linéaire de base. Enfin, nous fixons maxI le nombre d’itérations permettant de déterminer la convergence ou non de l’algorithme.

A noter qu’aussi bien pendant la phase d’initialisation que les autres étapes t, l’algorithme filtre les poids n’étant pas dans l’intervalle [B_{inf}, B_{sup}]. Ainsi, avant chaque déploiement de nouvelle étape, l’algorithme recherche les poids W ^t (égal à p s’il s’agit de la première itération) \notin [B_{inf}, B_{sup}] et les fixes aux bornes. Ainsi,

– si \exists i \in E, W_i ^t \mbox{ou } p_i < B_{inf} \Rightarrow W_i ^t \mbox{ou } p_i = B_{inf}

– si \exists i \in E, W_i ^t \mbox{ou } p_i > B_{sup} \Rightarrow W_i ^t \mbox{ou } p_i = B_{sup}

Donc, à l’itération t \leq maxI, nous avons tout d’abord la procédure de sélection des poids, seules ceux non réajustés en fonction des bornes sont conservés pour les calculs (la restriction sera indiquée par un astérisque). Puis la procédure de vérification visant à comparer les marges recalculées avec les pondérations W ^{t - 1} et les marges attendues M. Si elles ne sont pas toutes égales à \epsilon prés nous déployons la procédure ci-dessous pour l’étape en cours:

\begin{cases} \mbox{a - Calculer: } \lambda ^t = (\mathbf{X *} ' \cdot p * \cdot \mathbf{X *}) ^{-1} \cdot (p * \cdot \mathbf{X *} - M + W * ^{t - 1} \cdot p * \cdot \mathbf{X *}) \\ \mbox{b - Calculer: } W * ^t = 1 + \mathbf{X *} \cdot \lambda ^t \end{cases}

Les pondérations varient donc dans l’intervalle [B_{inf}, B_{sup}] puisque l’algorithme réajuste aux bornes les pondérations p puis W. Par conséquent la méthode, bien que linéaire à la base, ne génère plus de valeurs négatives.

Concernant la convergence, il est évident que le choix des bornes vont avoir une importance particulière puisque l’algorithme ne traite que les poids qui ne seront pas réajusté et donc qui appartiennent à [B_{inf}, B_{sup}]. Outre ce détail, la convergence de l’algorithme devrait être assurée étant donné qu’il se base sur la régression linéaire qui ne souffre pas de complication pour estimer les paramètres.

4) La méthode logistique

La méthode logistique est un algorithme itératif de résolution de \lambda.

La phase d’initialisation consiste à fixer les bornes B_{inf}, B_{sup} de l’intervalle de variation de W, puis à initialiser:

\begin{cases} c = \frac{B_{inf} - B_{sup}}{(1 - B_{inf}) \cdot (B_{sup} - 1)} \\ W=\frac{B_{inf} \cdot (B_{sup} - 1) + B_{sup} \cdot (1 - B_{inf})}{B_{sup} - B_{inf}} \\ \lambda = 0 \end{cases}

Enfin, nous fixons maxI le nombre d’itérations permettant de déterminer la convergence ou non de l’algorithme.

A noter lors de chaque nouvelle étape t, l’algorithme filtre les poids n’étant pas dans l’intervalle [B_{inf}, B_{sup}]. Avant chaque déploiement de nouvelle étape, l’algorithme recherche les poids W ^t \notin [B_{inf}, B_{sup}] et les fixes aux bornes, ainsi,

– si \exists i \in E, W_i ^t < B_{inf} \Rightarrow W_i ^t = B_{inf}

– si \exists i \in E, W_i ^t > B_{sup} \Rightarrow W_i ^t = B_{sup}

De plus, la méthode logistique modifie la table des marges M en lui soustrayant les totaux pondérés des observations retirées de \mathbf{X},

M * = M - (W * ^{t - 1} \cdot p *)' \cdot \mathbf{X *}

Donc, à l’itération t \leq maxI, nous avons tout d’abord la procédure de sélection des poids, seules ceux non réajustés en fonction des bornes sont conservés pour les calculs (la restriction sera indiquée par un astérisque). Puis la procédure de vérification visant à comparer les marges recalculées avec les pondérations W ^{t - 1} et les marges attendues M *. Si elles ne sont pas toutes égales à \epsilon prés nous déployons la procédure ci-dessous pour l’étape en cours:

\begin{cases} \mbox{a - Calculer: } \lambda ^t = \lambda ^{t - 1} - [(\mathbf{X *} \cdot p * \cdot W * ^{t - 1})' \cdot \mathbf{X *}] ^{-1} \cdot (\mathbf{X *} ' \cdot p * \cdot W * ^{t - 1} - M * + p * ^t \cdot \mathbf{X *}) \\ \mbox{b - Calculer: } U ^t = e ^{c \cdot \mathbf{X *} \cdot \lambda ^t} \\ \mbox{c - Calculer: } W ^t = \frac{B_{inf} \cdot (B_{sup} - 1) + B_{sup} \cdot (1 - B_{inf}) \cdot U ^t }{B_{sup} - 1 + (1 + B_{inf}) \cdot U ^t} \end{cases}

A l’instar de la méthode du Raking Ratio, l’intervention de l’exponentielle va permettre d’assurer la positivité des coefficients.

Concernant la convergence, ce qui fait la force de la méthode logistique va également être sa faiblesse. La fonction exponentielle augmente les chances de rapidement faire exploser les pondérations W, néanmoins les bornes [B_{inf}, B_{sup}] présentent l’avantage de fournir un intervalle permettant de diminuer les chances d’être dans cette situation.

\bullet Exemple:

Soit l’échantillon suivant E:

add

, obtenu suite à un sondage mené sur une sous-population de taille n = 15 et tirée aléatoirement selon une méthode quelconque depuis une population d’étude de taille N = 230. E se décrit donc au travers des variables X^1 qualitative et X ^2 quantitative et de la colonne p = pond qui présente les différents poids de sondage.

Le travail bibliographique mené sur la population d’étude nous permet d’avoir sous la main la table des marges réelles suivantes M:

add

Nous avons donc,

\mathbf{X} = \begin{pmatrix} 1 & 0 & 150 \\ 1 & 0 & 360 \\ 1 & 0 & 250 \\ 0 & 1 & 300 \\ 1 & 0 & 320 \\ 0 & 1 & 350 \\ 1 & 0 & 160 \\ 0 & 1 & 220 \\ 1 & 0 & 210 \\ 0 & 1 & 160 \\ 1 & 0 & 270 \\ 0 & 1 & 250 \\ 1 & 0 & 210 \\ 1 & 0 & 180 \\ 1 & 0 & 160 \end{pmatrix}

et, 

\mathbf{M} = \begin{pmatrix} 140 & 90 & 47000 \end{pmatrix}

p = (10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20, 20, 20, 20)

1) Redressement des pondérations par la méthode linéaire:

Pour l’application de la méthode linéaire il faut, dans un premier temps, calculer:

\lambda = [(\mathbf{X} \times p)' \times \mathbf{X}] ^{-1} \times (\mathbf{M} - p' \times \mathbf{X})

= \begin{pmatrix} 150 & 0 & 33000 \\ 0 & 80 & 19100 \\ 33000 & 19100 & 12663000 \end{pmatrix} ^{-1} \times \begin{pmatrix} - 10 \\ 10 \\ -5100 \end{pmatrix}

= \begin{pmatrix} 1.313431213 \\ 1.622719858 \\ -0.06273172 \end{pmatrix}

Les multiplicateurs de Lagrange \lambda étant déterminés, il nous reste plus qu’à redresser les p en W grâce à \lambda. Par conséquent,

W = 1 + \mathbf{X} \times \lambda = 1 + \begin{pmatrix} 0.37245586 \\ -0.944910772 \\ \vdots \\ 0.18426022 \\ 0.30972366 \end{pmatrix} = \begin{pmatrix} 1.37245539 \\ 0.05508923 \\ 0.74513817 \\ 0.74076820 \\ 0.30601612 \\ 0.42710960 \\ 1.30972366 \\ 1.24262198 \\ 0.99606506 \\ 1.61901231 \\ 0.61967472 \\ 1.05442681 \\ 0.99606506 \\ 1.18426022 \\ 1.30972366 \end{pmatrix}

2) La méthode Raking Ratio

Nous initialisons \lambda ^0 = (0, \cdots, 0), W ^0 = p et fixons maxI à un nombre quelconque pour cet exemple. Lançons la procédure itérative.

Pour l’itération N°1, nous avons donc:

\lambda ^1 = \lambda ^0 - [(\mathbf{X} \cdot W ^{t - 1})' \cdot \mathbf{X}] ^{-1} \cdot (\mathbf{X}' \cdot W ^{t - 1} - \mathbf{M})

= \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} - \begin{pmatrix} 150 & 0 & 33000 \\ 0 & 80 & 19100 \\ 33000 & 19100 & 12663000 \end{pmatrix} ^{-1} \times \begin{pmatrix} 10 \\ -10 \\ 5100 \end{pmatrix}

= \begin{pmatrix} 1.313431213 \\ 1.622719858 \\ -0.006273172 \end{pmatrix}

Et donc,

W ^1 = p \times e ^{\mathbf{X} \times \lambda ^1} = \begin{pmatrix} 10 \\ 10 \\ \vdots \\ 20 \\ 20 \end{pmatrix} \times \begin{pmatrix} 1.4512937 \\ 0.3887143 \\ \vdots \\ 1.2023287 \\ 1.3630484 \end{pmatrix} = \begin{pmatrix} 14.512937 \\ 3.887143 \\ \vdots \\ 24.046573 \\ 27.260968 \end{pmatrix}

L’itération N° 1 est achevée, nous regardons si les marges recalculées grâce à W ^{1} s’approche de celles de \mathbf{M},

\mathbf{X} ' \times W ^1 - \mathbf{M} = \begin{pmatrix} 9.59985 >>>> 0 \\ 7.10767 >>>> 0 \\ 4161.32358 >>>> 0 \end{pmatrix}

Pour cette première itération la distance entre les marges recalculées et les marges attendues restent considérables. Par conséquent, nous poursuivons la procédure itérative jusqu’à ce nous ayons convergé vers les marges attendues.

Pour l’itération N° 2 nous obtenons la différence entre marges recalculées avec W ^2 et marges attendues M suivante,

\mathbf{X} ' \times W ^2 - M = \begin{pmatrix} 0.4819199 > 0 \\ 0.3975174 > 0 \\ 245.1318572 >>>> 0 \end{pmatrix}

Pour l’itération N° 3 nous obtenons la différence entre marges recalculées avec W ^3 et marges attendues M suivante,

\mathbf{X} ' \times W ^3 - M = \begin{pmatrix} 0.002740598 > 0 \\ 0.002427357 > 0 \\ 1.501063415 > 0 \end{pmatrix}

Pour l’itération N° 4 nous obtenons la différence entre marges recalculées avec W ^4 et marges attendues M suivante,

\mathbf{X} ' \times W ^4 - M = \begin{pmatrix} 0.000000171574 \approx 0 \\ 0.0000001062141 \approx 0 \\ 0.00006496802 \approx 0 \end{pmatrix}

L’algorithme a donc convergé en quatre itérations, nous obtenons donc:

W = W ^4 = \begin{pmatrix} 14.366193 \\ 3.031181 \\ 6.848029 \\ 6.499062 \\ 4.076827 \\ 4.487068 \\ 13.340253 \\ 23.512650 \\ 18.420695 \\ 36.674757 \\ 11.809740 \\ 18.826463 \\ 18.420695 \\ 23.005880 \\ 26.680507 \end{pmatrix}

3) La méthode linéraire tronquée

Pour la phase d’initialisation nous fixons B_{inf} = 0.3, B_{sup} = 1.8 pour la variation des W. Nous reprenons ensuite les pondérations calculées par la méthode linéaire de base et initialisons,

W ^1= (1.37245539, \cdots, 1.30972366)

Enfin, nous fixons maxI à un nombre quelconque pour cet exemple. Lançons la procédure itérative.

Nous démarrons l’itération N° 1. L’étape préliminaire de chaque itération demande de prendre W et de borner ses valeurs en fonction de B_{inf}, B_{sup}. Par conséquent nous remarquons que seul: W_2 = 0.05508923 est hors des bornes, nous fixons donc W_2 = B_{inf} = 0.3. Pour la suite de l’application, nous poserons W*, \mathbf{X} *, p*, respectivement, le vecteur des poids redressés, la matrice de données et les pondérations initiales privés de l’observation N° 2.

Nous avons alors,

\lambda ^1 = (\mathbf{X *} ' \cdot p * \cdot \mathbf{X *}) ^{-1} \cdot (p * \cdot \mathbf{X *} - M + W * ^{t - 1} \cdot p * \cdot \mathbf{X *})

= \begin{pmatrix} 140 & 0 & 29400 \\ 0 & 80 & 19100 \\ 29400 & 19100 11367000 \end{pmatrix} ^{-1} \times \begin{pmatrix} -3 \\ 10 \\ -2580 \end{pmatrix}

= \begin{pmatrix} 1.417836686 \\ 1.761307525 \\ -0.006853644 \end{pmatrix}

Nous pouvons désormais calculer les nouvelles pondérations,

W * ^1 = 1 + \mathbf{X *} \times \lambda ^1 = 1 + \begin{pmatrix} 0.38979007 \\ -0.29557433 \\ \vdots \\ 0.18418075 \\ 0.32125363 \end{pmatrix} = \begin{pmatrix} 1.3897901 \\ 0.7044257 \\ \vdots \\ 1.1841808 \\ 1.3212536 \end{pmatrix}

Enfin,

W ^1 = (W_1 * ^1, W_2 ^1 = B_{inf}, W_3 * ^1, \cdots, W_15 * ^1) = (1.3897901, 0.3, 0.7044257, 0.7052143, 0.2246706, \cdots, 1.3212536)

Nous constatons que W_5 ^1 = 0.2246706 < B_{inf}, par conséquent nous allons recourir à une nouvelle itération qui visera à borner W_5 ^1 tout en cherchant l’égalité des marges observées et attendues.

Pour l’itération N° 2, nous avons donc deux pondérations qui sont hors des bornes B_{inf}, B_{sup} donc nous ne traitons que les observations N° 1, 3, 4, 6, 7, \cdots, 15. Les observations N° 2, 6 étant fixées à B_{inf}. Nous obtenons

W ^2 = (1.393467 ,0.3, 0.6900262, 0.6943387, 0.3, \cdots, 1.3228347)

, dont tous les poids sont bien dans l’intervalle [B_{inf}, B_{sup}]. Enfin, pour savoir si nous devons procéder à une nouvelle itération, nous vérifions si les marges recalculées à partir de W ^2 sont suffisamment proches de celles attendues dans M,

\mathbf{X} ^t \times W ^2 - M = \begin{pmatrix} 140 \\ 90 \\ 47000 \end{pmatrix} - \begin{pmatrix} 140 \\ 90 \\ 47000 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}

Par conséquent, l’algorithme a bien convergé en deux itérations et sous contrainte de respecter l’intervalle fixé pour les pondérations possibles de W. Nous avons donc,

W = W ^2 = \begin{pmatrix} 1.3931467 \\ 0.3 \\ 0.6900262 \\ 0.6943387 \\ 0.3 \\ 0.3427784 \\ 1.3228347 \\ 1.2568351 \\ 0.9712744 \\ 1.6787074 \\ 0.5494021 \\ 1.0458989 \\ 0.9712744 \\ 1.1822106 \\ 1.3228347 \end{pmatrix}

4) La méthode logistique

Pour la phase d’initialisation nous fixons B_{inf} = 0.3, B_{sup} = 1.8 pour la variation des W. Ensuite, nous fixons maxI à un nombre quelconque pour cet exemple. Nous devons également calculer,

\begin{cases} c = \frac{0.3 - 1.8}{(1 - 0.3) \cdot (1.8 - 1)} = 2.678571 \\ W ^0 = \frac{0.3 \cdot (1.8 - 1) + 1.8 \cdot (1 - 0.3)}{1.8 - 0.3} = 1 \\ \lambda = 0 \end{cases}

Enfin, nous vérifions si les poids W intialisés ne suffisent pas à la convergence de l’algorithme,

\mathbf{X} ' \times W ^0 - \mathbf{M} = \begin{pmatrix} 150 \\ 80 \\ 52100 \end{pmatrix} - \begin{pmatrix} 140 \\ 90 \\ 47000 \end{pmatrix} = \begin{pmatrix} 10 >>>> 0 \\ -10 <<<< 0 \\ 5100 >>>> 0 \end{pmatrix}

Ce n’est pas le cas, nous lançons la procédure itérative.

Nous démarrons l’itération N° 1. L’étape préliminaire de chaque itération demande de prendre W et de borner ses valeurs en fonction de B_{inf}, B_{sup}. Pour cette itération, les poids W valent tous égaux à 1 et sont donc inclus dans l’intervalle. La table des marges M est modifiée dans le cas où au moins un poids W n’appartient pas à l’intervalle fixé. Ce qui n’est pas le cas ici, par conséquent,

M * = M

Et \mathbf{X *} = \mathbf{X}, W * = W et p * = p pour la même raison.

Nous pouvons donc calculer,

\lambda ^1 = \lambda ^0 - [(\mathbf{X *} \cdot p * \cdot W * ^0)' \cdot \mathbf{X *}] ^{-1} \cdot (\mathbf{X *} ' \cdot p * \cdot W * ^0 - M * + p * \cdot \mathbf{X *})

= \begin{pmatrix} 0 \\ 0 \\ \end{pmatrix} - \begin{pmatrix} 150 & 0 & 33000 \\ 0 & 80 & 19100 \\ 3300 & 19100 & 12663000 \end{pmatrix} ^{-1} \times \begin{pmatrix} 10 \\ -10 \\ 5100 \end{pmatrix}

= \begin{pmatrix} 2.7612571 \\ 3.42642009 \\ -0.01330438 \end{pmatrix}

Nous pouvons ainsi déterminer,

U ^1 = e ^{c \times \mathbf{X *} \times \lambda ^1} = e ^{2.678571 \times \begin{pmatrix} 0.76560211 \\ -2.02831713 \\ \vdots \\ 0.36647079 \\ 0.63255834 \end{pmatrix}} = \begin{pmatrix} 7.77349559 \\ 0.00437000 \\ \vdots \\ 2.66877134 \\ 5.44312147 \end{pmatrix}

Et,

W ^1 = \frac{0.3 \times (1.8 - 1) + 1.8 \times (1 - 0.3) \times \begin{pmatrix} 7.77349559 \\ 0.00437000 \\ \vdots \\ 2.66877134 \\ 5.44312147 \end{pmatrix}}{(1.8 - 1 + (1 - 0.3) \times \begin{pmatrix} 7.77349559 \\ 0.0043700 \\ \vdots \\ 2.66877134 \\ 5.44312147 \end{pmatrix}} = \begin{pmatrix} 1.3552801 \\ 0.3976475 \\ \vdots \\ 1.1835689 \\ 1.3009790 \end{pmatrix}

Nous constatons que l’ensemble des poids W ^1 sont dans l’intervalle [B_{inf}, B_{sup}]. Il nous reste à vérifier si les marges recalculées à partir de W ^1 sont suffisamment proches de celle attendues M. Nous avons,

\mathbf{X} ' \times W ^1 - M = \begin{pmatrix} 145.67357 \\ 89.45847 \\ 49117.01517 \end{pmatrix} - \begin{pmatrix} 140 \\ 90 \\ 47000 \end{pmatrix} = \begin{pmatrix} 5.6735740 >>>> 0 \\ -0.5415341 < 0 \\ 2117.0151703 >>>> 0 \end{pmatrix}

Etant donné la trop grande différence entre les deux marges, nous poursuivons l’algorithme itératif jusqu’à ce que les deux conditions soient remplies.

Nous obtenons donc:

– Pour l’itération N°2, W ^2 = (1.3959632, 0.3438108 \cdots, 1.1816716, 1.3296580), nous avons zéro poids hors de l’intervalle [0.3, 1.8] et la différence entre marges recalculées à partir de W ^2 et marges attendues M donne,

\mathbf{X} ' \times W ^2 - M = \begin{pmatrix} 140.38621 \\ 87.94418 \\ 46971.57196 \end{pmatrix} - \begin{pmatrix} 140 \\ 90 \\ 47000 \end{pmatrix} = \begin{pmatrix} 0.3862149 > 0 \\ -2.0558175 <<<< 0 \\ -28.4280384 <<<< 0 \end{pmatrix}

– Pour l’itération N°3, W ^3 = (1.4160877, 0.3338588 \cdots, 1.1895082, 1.3467227), nous avons zéro poids hors de l’intervalle [0.3, 1.8] et la différence entre marges recalculées à partir de W ^3 et marges attendues M donne,

\mathbf{X} ' \times W ^3 - M = \begin{pmatrix} 139.8655 \\ 88.5988 \\ 46823.8668 \end{pmatrix} - \begin{pmatrix} 140 \\ 90 \\ 47000 \end{pmatrix} = \begin{pmatrix} -0.1345485 < 0 \\ -1.4012000 < 0 \\ -176.1332391 < 0 \end{pmatrix}

\vdots

– Pour l’itération N°22, W ^{22} = (1.4398983, 0.3263306, \cdots, 1.2026114, 1.3682038), nous avons zéro poids hors de l’intervalle [0.3, 1.8] et la différence entre marges recalculées à partir de W ^{22} et marges attendues M donne,

\mathbf{X} ' \times W ^{22} - M = \begin{pmatrix} 139.99999 \\ 89.99994 \\ 46999.99046 \end{pmatrix} - \begin{pmatrix} 140 \\ 90 \\ 47000 \end{pmatrix} = \begin{pmatrix} -0.000009395427 \approx 0 \\ -0.00006018689 \approx 0 \\ -0.000953811 \approx 0 \end{pmatrix}

Par conséquent, l’algorithme converge au bout de 22 itérations et les poids redressés sont,

W = W ^{22} = \begin{pmatrix} 1.4398983 \\ 0.3263306 \\ 0.61798952 \\ 0.6268710 \\ 0.3685640 \\ 0.4126877 \\ 1.3682038 \\ 1.3004038 \\ 0.9285075 \\ 1.6467703 \\ 0.5116785 \\ 1.03330436 \\ 0.9285075 \\ 1.2026114 \\ 1.3682038 \end{pmatrix}

\bullet Application informatique:

– sous SAS: http://www.insee.fr/fr/methodes/default.asp?page=outils/calmar/accueil_calmar.htm

– sous R: http://maths.cnam.fr/IMG/txt/Exemple_calage_avec_R_cle4611bd.txt

\bullet Bibliographie:

– Les techniques de sondages de Pascal Ardilly

– Probabilité, analyse des données et statistique de Gilbert Saporta

– La macro CALMAR. Redressement d’un échantillon par calage sur Marges d’ Olivier Sautory

Les plans d’expérience

add\bullet Introduction:

Une expérience est une opération menée sous des hypothèses fixées et entièrement contrôlées dans le but de découvrir un effet inconnu, de tester une relation ou encore de démontrer une loi connue et donc clarifier la relation entre ces hypothèses et le résultat de l’expérience. L’experience devient une expérience aléatoire si les résultats possibles peuvent être décrit au travers d’une probabilité de réalisation bien définie. Dans le cadre de la présence de facteurs, nous parlerons d’expérience factorielle. Le plan d’expériences est alors construit en fonction de l’objectif de l’étude, des facteurs intervenant et du coût d’opération.

Les plans d’expérience représentent donc la seconde approche de collecte des données avec les plans de sondage. Ils ont pour objectif de proposer un modèle reliant les variables explicatives \mathbf{X} (appelés facteurs dans le jargon des plans d’expériences) à la variable à expliquer \mathbf{Y} et visant ainsi à optimiser l’étude de l’intéraction entre \mathbf{X} avec \mathbf{Y}. Le grand intérêt des plans d’expérience étant notamment de déterminer le nombre d’expériences minimum et optimal à effectuer afin de concilier coût de calculs et fiabilité des résultats.

Chaque plan d’expérience possède ses hypothèses d’utilisation et sa méthodologie à suivre.

\bullet Le problème:

L’objectif est d’estimer au mieux le vecteur \beta de paramètres du modèle linéaire:

\mathbf{Y} = \mathbf{X} \cdot \mathbf{\beta} + \mathbf{\epsilon}

, ce qui revient à obtenir des estimateurs de variance minimale et dépendant exclusivement de \mathbf{X}. Pour se faire, deux principaux critères sont utilisés:

– Le critère du déterminant maximal ou D-optimalité: max | \mathbf{X' X} |

– Le critère de la somme des variances minimale ou A-optimalité: min (Tr (\mathbf{X ' X} ^{-1}))

Concernant les facteurs qualitatifs, il est intéressant de noter qu’une condition d’orthogonalité nécessaire et suffisante pour les modèles à effets principaux sans interaction est que le plan soit équlibré ce qui revient à dire que nous sommes en présence du même nombre de mesures pour \mathbf{Y} quelque soit le croisement de nos facteurs.

Remarquons que pour des facteurs ordinaux ce problème revient à étudier la corrélation deux à deux des différents facteurs. Dans le cas de non-corrélation nous avons orthogonalité entre les facteurs ce qui implique un plan optimal. Ce fait n’est pas toujours respecté, nous nous intéressons alors à l’isovariance par rotation soit l’étude de la variance de la prédiction de \mathbf{Y} en un point de \mathbf{X}.

Une fois le plan d’expérience définit, il suffit alors de mesurer \mathbf{Y} pour les diférentes combinaisons de facteurs retenues et appliquer une analyse statistique (en générale l’analyse de variance) afin d’obtenir des résultats fiables. En général, les mesures sont prises une fois par combinaison de facteurs retenus, néanmoins il est tout à fait possible de prendre plusieurs mesures si l’expérimentateur juge plus important de diminuer le nombre d’expériences à effectuer au détriment du nombre de mesures à prendre.

\bullet Les différents plans:

1) Les plans carrés latins

Histoire: Première trace d’utilisation en 1624.

Principe: Les plans carrés latins sont des plans orthogonaux qui consistent à considérer trois facteurs qualitatifs (X ^1, X ^2, X ^3) à K \geq 2 groupes et de construire la matrice du plan d’expérience de dimension K \times K présentant les différentes combinaisons des trois facteurs croisés.

Critique: Les plans carrés latins permettent de réduire le nombre d’expériences de K ^3 à K ^2. Ils ne peuvent estimer que les effets principaux. Ce sont des plans orthogonaux D-optimal.

Exemple: Soit les trois facteurs (X ^1, X ^2, X ^3) à trois groupes tel que X ^1 \in \lbrace A1, B1, C1 \rbrace, X ^2 \in \lbrace A2, B2, C2 \rbrace et X ^3 \in \lbrace A3, B3, C3 \rbrace. La matrice du plan d’expérience est alors de dimensions 3 \times 3 et vaut:

\begin{pmatrix} X ^1 = / X ^2 = & A2 & B2 & C2 \\ A1 & A3 & B3 & C3 \\ B1 & B3 & C3 & A3 \\ C1 & C3 & A3 & B3 \end{pmatrix}

, la première ligne et la première colonne correspondent respectivement aux modalités prises par les facteurs X ^2, X ^1 et les cellules restantes la valeur du facteur X ^3 quand X ^1 = k_1, X ^2 = k_2.

2) Les plans carrés gréco-latins

Histoire: Première trace d’utilisation en 1782 par Leonhard Euler.

Principe: Les plans carrés gréco-latins sont des plans orthogonaux qui consistent à considérer quatre facteurs qualitatifs (X ^1, X ^2, X ^3, X ^4) à K \geq 2 groupes en superposant les deux plans carrés latins sur (X ^1, X ^2, X ^3) et (X ^1, X ^2, X ^4).

Critique: Les plans gréco-latins permettent de réduire le nombre d’expériences de K ^4 à K ^2. Il existe des plans carrés gréco-latins pour toutes dimensions à l’exception de K = 1, 2 et 6. Ce sont des plans orthogonaux D-optimal.

Exemple: Soit les quatre facteurs (X ^1, X ^2, X ^3, X ^4) à K = 3 tel que X ^1 \in \lbrace A1, B1, C1 \rbrace, X ^2 \in \lbrace A2, B2, C2 \rbrace, X ^3 \in \lbrace A3, B3, C3 \rbrace et X ^4 \in \lbrace A4, B4, C4 \rbrace. La matrice du plan d’expérience est alors,

\begin{pmatrix} X ^1 = / X ^2 = & A2 & B2 & B3 \\ A1 & A3 A4 & B3 B4 & C3 C4 \\ B1 & B3 C4 & C3 A4 & A3 B4 \\ C1 & C3 B4 & A3 C4 & B3 A4 \end{pmatrix}

, la première ligne et la première colonne correspondent respectivement aux modalités prises par les facteurs X ^2, X ^3 et les cellules restantes la valeur du croisement des facteurs X ^3, X ^4 quand X ^1 = k_1, X ^2 = k_2.

3) Les plans complètement randomisés

Histoire: Première trace d’utilisation en 1920 par Ronald Aylmer Fisher.

Principe: Les plans complètement randomisés consistent à considérer un ou plusieurs facteurs qualitatifs (X ^1, \cdots, X ^P) à K_1, \cdots, K_P \geq 2 groupes respectifs et attribuer aléatoirement pour chaque croisement de groupes des différents facteurs un nombre n_{k_1, \cdots, n_{k_P}} de mesures à prendre au sein de l’ensemble des unités statistiques supposées uniformes.

Si toutes les combinaisons apparaissent au moins une fois alors le plan est considéré complet, sinon incomplet.

Critique: Globalement il s’agit du plan d’expérience se rapprochant le plus d’un tirage aléatoire simple. Pas forcément le plus économique, néanmoins il permet de considérer des facteurs à nombre de modalités différents.

Exemple: Soit les deux facteurs (X ^1, X ^2) à, respectivement, K_1 = 3 et K_2 = 2 groupes. Nous décidons de tirer aléatoirement un effectif pour les différentes combinaisons de groupes des facteurs et obtenons n_{k_1 = 1, k_2 = 1} = 2, n_{1,2} = 3, n_{2,1} = 1, n_{2,2} = 5, n_{3,1} = 3, n_{3,2} = 4 d’observations. La matrice du plan d’expérience est alors,

\begin{pmatrix} X ^1 & X ^2 & Effectif \\ g_1 & g_1 \\ \vdots & & \times 2 \\ g_1 & g_2 \\ \vdots & & \times 3 \\ g_2 & g_1 \\ \vdots & & \times 1 \\ g_2 & g_2 \\ \vdots & & \times 5 \\ g_3 & g_1 \\ \vdots & & \times 3 \\ g_3 & g_2 \\ \vdots & & \times 4 \end{pmatrix}

4) Les plans randomisés à blocs

Histoire: Première trace d’utilisation en 1920 par Ronald Aylmer Fisher.

Principe: Les plans randomisés à blocs consistent à considérer des combinaisons de groupes de un ou plusieurs facteurs qualitatifs (X ^1, \cdots, X ^P) à K_1, \cdots, K_P \geq 2 groupes respectifs. Ces q = K_1 \times \cdots \times K_P combinaisons donnent q! permutations possibles et sont alors arrangées et tirées r (fixé au prélable) fois aléatoirement avec remise pour constituer les blocs du plan d’expérience.

Critique: Les plans randomisés par blocs permettent de réduire le nombre d’expériences de q! à r. De plus, ils permettent de considérer des facteurs à nombre de modalités différent.

Exemple: Soit les deux facteurs (X ^1, X ^2) à, respectivement, K_1 = 3 et K_2 = 2 groupes. Nous avons alors les combinaisons suivantes A = (X ^1 = g_1, X ^2 = g_1), B = (X ^1 = g_1, X ^2 = g_2), C = (X ^1 = g_2, X ^2 = g_1), D = (X ^1 = g_2, X ^2 = g_2), E = (X ^1 = g_3, X ^2 = g_1), F = (X ^1 = g_3, X ^2 = g_2) ce qui constitue 3 \times 2 = 6 combinaisons possibles et donc 6! = 720 permutations. La matrice du plan d’expérience est alors, pour r = 6 blocs construient par tirage aléatoire avec remise sur les différentes permutations de combinaisons,

\begin{pmatrix} \mbox{bloc 1 = } & E & D & C & B & A & F \\ \mbox{bloc 2 = } & D & C & B & A & E & F \\ \mbox{bloc 3 = } & E & B & A & C & D & F \\ \mbox{bloc 4 = } & E & D & C & B & A & F \\ \mbox{bloc 5 = } & A & B & C & F & E & D \\ \mbox{bloc 6 = } & E & C & D & A & B & F \end{pmatrix}

5) Les plans factoriels complets

Histoire: ND.

Principe: Les plans factoriels complets sont des plans orthogonaux qui consistent à considérer un ou plusieurs facteurs (X ^1, \cdots, X ^P) ordinaux en posant, par convention, -1 la valeur minimale pour le facteur X ^p et +1 pour la valeur maximale. Le plan d’expériences se construit à partir des 2 ^P expériences possibles soient les combinaisons des couples (-1, 1) des différents facteurs considérés. Si nous voulons intégrer les intéractions, il suffit de multiplier les colonnes des facteurs seuls entre eux pour construire celles des intéractions.

Si les combinaisons de niveaux sont présentes le même nombre de fois alors le plan factoriel complet est dit équilibré.

Les plans fractionnaires de type 2 ^{P - k} sont dérivés des plans factoriels complets et proposent une solution au problème suivant: proposer un plan d’expériences permettant d’étudier P facteurs à partir d’un nombre d’expériences moindres. L’approche est alors de considérer le plan factoriel complet associé au nombre maximum de facteurs compatibles avec le nombre d’expériences. Ce sous-ensemble, épuré de k facteurs et qui constituent les facteurs de « base », permet alors de construire les nouveaux facteurs qui sont le fruit du produit entre les facteurs de « base ». Nous dirons qu’ils sont confondus avec l’intéraction considérée. Ce type de plan se nomme également plan de Box et Hunter ou de Plackett et Burman (si le nombre d’expérience n’est plus un multiple de 2 mais de 4) et constituent les plans de criblage principalement adaptés aux études de faisabilité.

Critique: Ce sont des plans orthogonaux D-optimal et A-optimal.

Exemple: Soit les trois facteurs (X ^1, X ^2, X ^3) à deux niveaux chacun. Nous avons alors 2 ^3 = 8 expériences possibles. La matrice du plan d’expérience factoriel complet pour les facteurs d’ordre 1 et le facteur d’ordre 2 X ^1 \times X ^2 est alors,

\begin{pmatrix} Experience & X ^1 & X ^2 & X ^3 & X ^1 \times X ^2 \\ 1 & -1 & -1 & -1 & +1 \\ 2 & +1 & -1 & -1 & -1 \\ 3 & -1 & +1 & -1 & -1 \\ 4 & +1 & +1 & -1 & +1 \\ 5 & -1 & -1 & +1 & +1 \\ 6 & +1 & -1 & +1 & -1 \\ 7 & -1 & +1 & +1 & -1 \\ 8 & +1 & +1 & +1 & +1 \end{pmatrix}

6) Les plans composites

Histoire: ND.

Principe: Les plans composites sont des plans qui ne sont pas orthogonaux et qui consistent à considérer un ou plusieurs facteurs ordinaux (X ^1, \cdots, X ^P) seuls, à 3 niveaux, avec leur intéraction d’ordre 2 (soit le croisement deux à deux des facteurs considérés). Il faut poser, par convention, -1 la valeur minimale pour le facteur X ^p0 sa valeur médiane et +1 pour la valeur maximale. Les plans composites se placent dans une approche géométrique du problème en cherchant à visualiser le carré (pour deux facteurs), le cube (pour trois facteurs) ou l’hyper-cube (pour plus de trois facteurs) construit à partir du codage ci-dessus. Les P facteurs sont placés aux extrémités et les intéractions, d’ordre 2, le sont à une distance \alpha (fixée) des extrémités. Le plan d’expérience se construit à l’instar d’un plan factoriel complet (ou fractionnaire) et considère pour les facteurs d’ordre 2 la matrice d’experience de taille [2 \cdot (P - k)] \times (P - k) (où on soustrait à P le nombre k de facteurs non considérés comme de « base » si nous sommes en présence d’un plan fractionnaire, k = 0 pour le plan factoriel complet) de forme:

\begin{pmatrix} Experience & X ^1 & X ^2 & X ^3 & \cdots & X ^{P - k - 1} & X ^{P - k} \\ 2 ^{P - k} + 1 & -\alpha & 0 & 0 & \cdots & 0 & 0 \\ 2 ^{P - k} + 2 & +\alpha & 0 & 0 & \cdots & 0 & 0 \\ 2 ^{P - k} + 3 & 0 & -\alpha & 0 & \cdots & 0 & 0 \\ 2 ^{P - k} + 4 & 0 & +\alpha & 0 & \cdots & 0 & 0 \\ \vdots \\ 2 ^{P - k} + 2 \cdot (P - k - 1) & 0 & 0 & 0 & \cdots & 0 & -\alpha \\ 2 ^{P - k} + 2 \cdot (P - k) & 0 & 0 & 0 & \cdots & 0 & +\alpha \end{pmatrix}

, également complété par L0 lignes nulles représentant les points à l’origine de la figure géométrique.

Les principales caractéristiques des plans composites étant l’isovariance par rotation et la possibilité d’expérimentation séquentielle qui revient à augmenter un plan factoriel fractionnaire de criblage en lui ajoutant des points au centre et d’autres points pour estimer les autres effets. Les plans composites appartiennent à la famille des plans pour surfaces de réponse du second degré. Le plan de Box-Behnken (de la même famille) représente une alternative au plan composite qui n’est pas isovariant par rotation mais adapté aux expériences difficiles à réaliser. Le plan d’expérience devient alors l’ensemble des permutations du triplet (-1, 0, +1) tel qu’au moins un facteur soit égal à 0 par expérience. Le plan de Doelhert, dont la stratégie repose sur la visualisation géométrique d’un polygone et plus d’un cube, peut également être cité comme l’un des plus connus.

Critique: Les plans composites ne sont pas orthogonaux mais permettent de travailler sur des facteurs à K = 3 niveaux.

Exemple: Soit les trois facteurs (X ^1, X ^2, X ^3) à trois niveaux chacun. Nous avons alors 2 ^3 = 8 expériences possibles pour le plan factoriel complet. La matrice du plan d’expérience se divise alors en 3 blocs, pour \alpha = 1, L0 = 2,

\begin{pmatrix} Experience & X ^1 & X ^2 & X ^3 \\ 1 & -1 & -1 & -1 \\ 2 & +1 & -1 & -1 \\ 3 & -1 & +1 & -1 \\ 4 & +1 & +1 & -1 \\ 5 & -1 & -1 & +1 \\ 6 & +1 & -1 & +1 \\ 7 & -1 & +1 & +1 \\ 8 & +1 & +1 & +1 \\ & & & \\ 9 & -1 & 0 & 0 \\ 10 & +1 & 0 & 0 \\ 11 & 0 & -1 & 0 \\ 12 & 0 & +1 & 0 \\ 13 & 0 & 0 & -1 \\ 14 & 0 & 0 & +1 \\ & & & \\ 15 & 0 & 0 & 0 \\ 16 & 0 & 0 & 0 \end{pmatrix}

7) Les plans asymétriques

Histoire: ND.

Principe: Les plans asymétriques consistent à considérer un ou plusieurs facteurs qualitatifs (X ^1, \cdots, X ^P) à, respectivement, K_1, \cdots, K_P \geq 2 groupes. Il existe plusieurs méthodes sous contrainte que le nombre d’expériences minimum doit être égal au nombre de paramètres à estimer  > \sum_{p = 1} ^P (K_p - 1) + 1. La plus populaire est celle de fusion qui consiste à considérer le plan fractionnaire pour P - k facteurs de « base » et à construire par combinaison de facteurs ceux dont nous avons besoin pour respecter le nombre de groupes respectifs voulus.

Critique: L’orthogonalité D-optimal est obtenu pour les plans asymétriques si le nombre d’expériences est un multiple des K_{p1}, K_{p2}, \forall p_1, p_2, ce qui est rare sauf dans le cas d’un plan complet.

Exemple: Nous voulons le plan pour trois facteurs, l’un de quatre niveaux et les deux autres à deux niveaux. Par la méthode de fusion. Nous partons du plan fractionnaire 2 ^{4 -1} = 8 suivant:

\begin{pmatrix} Experience & A & B & C & D \\ 1 & -1 & -1 & -1 & -1 \\ 2 & +1 & -1 & -1 & +1 \\ 3 & -1 & +1 & -1 & +1 \\ 4 & +1 & +1 & -1 & -1 \\ 5 & -1 & -1 & +1 & +1 \\ 6 & +1 & -1 & +1 & -1 \\ 7 & -1 & +1 & +1 & -1 \\ 8 & +1 & +1 & +1 & +1 \end{pmatrix}

Afin de construire notre facteur à quatre niveaux nous posons (C = -1, D = +1) \Rightarrow E = 1, (C = -1, D = +1) \Rightarrow E = 2, (C = +1, D = -1) \Rightarrow E = 3 et (C = +1, D = +1) \Rightarrow E = 1. Par conséquent nous obtenons le plan d’expérience suivant:

\begin{pmatrix} Experience & A & B & E \\ 1 & -1 & -1 & 1 \\ 2 & +1 & -1 & 2 \\ 3 & -1 & +1 & 2 \\ 4 & +1 & +1 & 1 \\ 5 & -1 & -1 & 4 \\ 6 & +1 & -1 & 3 \\ 7 & -1 & +1 & 3 \\ 8 & +1 & +1 & 4 \end{pmatrix}

\bullet Bibliographie:

– Statistique, dictionnaire encyclopédique de Yadolah Dodge

– Probabilité, analyse des données et Statistique de Gilbert Saporta

– Plan d’expérience pour l’estimation de la surface réponse de B. Govaerts

– Le document: http://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-m-modmixt4-plans.pdf

– Le document: Plans factoriels complets, plans fractionnaires. Cas des facteurs ayant deux modalités de Frédéric Bertrand

La théorie des sondages

add\bullet Introduction et justification:

La théorie des sondages a été mise en place afin d’offrir au statisticien (nous nous restreindrons au terme de sondeur, évidemment inclut dans tout statisticien maîtrisant intégralement les différentes facettes de sa discipline) des outils lui permettant de collecter l’information de manière optimale. Il existe deux types de techniques:

– le recensement,

– le plan de sondage.

Le premier consiste à déterminer une population d’étude et recueillir de manière exhaustive l’information. Le second consiste à tirer une sous-population de la population d’étude et recueillir l’information sur cette sous-population tout en garantissant qu’elle soit représentative de la population d’étude.

Il est évident que chaque sondeur rêve de pouvoir collecter exhaustivement l’information dont il a besoin, malheureusement cela entraîne deux principaux problèmes de taille:

– le temps de collecte,

– le coût pour collecter l’information.

Ces deux problèmes sont bien souvent liés et omniprésent pour toute étude de taille qui se respecte (Il reste évident que pour des études ciblant une région à la population très faible ou encore une population ultra-spécifique, le recensement reste de rigueur et n’est pas affligé par ces deux problématiques spécifiques). Afin d’y pallier, le sondeur s’orientera vers des méthodes moins contraigneuses telles que les plans de sondage. Les contraintes de temps et de coûts laissent alors place à une unique contrainte: la représentativité de la sous-population vis-à-vis de la population complète d’étude. Elle consiste à faire en sorte qu’une fois l’information collectée, les mesures calculées et associées (les estimateurs pour simplifier le texte) à l’information sur cette sous-population soient quasiment les même que celles que nous aurions obtenu sur la population totale à \epsilon, le plus faible possible, prés. C’est là tout l’enjeu des plans de sondages qui se déclinent au travers de six principales familles qui doivent être confrontées en fonction de la problématique et des informations à considérer (nommées régulièrement informations auxiliaires) dans l’objectif d’opter pour la méthode optimisant la représentativité de la sous-population.

Le problème de représentativité peut alors se décliner en trois problèmes plus accessibles à maîtriser:

– la sélection de l’échantillon d’individus (nommés unité de sondage ou unité statistique ou encore unité d’échantillonnage),

– l’agrégation de l’information, qui revient à une étape de détermination des estimateurs,

– la qualité de la procédure d’agrégation, soit la confiance que nous pouvons accorder aux résultats (nommé précision, terme à ne pas confondre avec le taux de bonne classification globale calculé de l’élaboration de modèles prédictifs en analyse supervisée et également nommé ainsi).

Ces trois problèmes sont fondamentalement interdépendant, impliquant que le sondage ne doit jamais être conçu autrement qu’en examinant parallèlement les trois points suivants:

– estimateur et précision dépendent fondamentalement de la méthode de sondage,

– la précision dépend de l’estimateur choisi,

– la méthode de sondage reste indépendante mais doit être choisie en fonction de la précision souhaitée et de la complexité de calcul des estimateurs d’inérêt.

En outre, déterminer la méthode de sélection de l’échantillon et la formulation de l’estimateur équivaut à déterminer la définition du plan de sondage. De plus, contrairement aux idées reçus, un biais nul n’est pas un critère infaillible, la priorité réside dans l’erreur quadratique moyenne.

Au final, l’enjeu revient à mettre en place un sondage qui va minimiser l’erreur totale soit la somme de l’erreur d’observation (que nous pouvons reformuler par le fait que la valeur mesurée pour une observation soit différente de la valeur réelle), des erreurs dues à une base de sondage incomplète ou comportant trop de non-réponses et de l’erreur d’échantillonnage (biais associé aux estimateurs et liés à la qualité du plan de sondage).

Cet enjeu est de taille, puisqu’au final nous ne connaissons jamais avec certitude la valeur des estimateurs d’intérêt. Le sondeur reste en aveugle et peut facilement être orienté vers de fausses pistes sans avoir conscience si son plan de sondage est mal composé. C’est bien pourquoi il faut prendre cette étape avec le plus de sérieux et de rigueur possible car elle reste l’étape de non-retour d’un projet d’étude.

\bullet Calcul du nombre de sujets nécessaires:

La question qu’il faut se poser avant de déterminer le plan de sondage conciliant précision et coût est la taille d’échantillon à enquêter permettant de coller aux objectifs fixés. Evidemment la théorie est là, les formules associées le sont également mais bien souvent le résultat du calcul du nombre de sujets nécessaires (NNS ou NSN) s’heurte à la réalité économique du projet.

Néanmoins, si nous restons dans un cadre purement théorique et en posant \epsilon le niveau de confiance fixé, m la marge d’erreur et \sigma ^2 la variance, les formules existantes du calcul du NSN en fonction des différents paramètres que nous souhaitons mesurer sont:

Pour l’estimation d’une proportion ou d’un pourcentage P:

n = \frac{\epsilon ^2 \cdot P \cdot (1 - P)}{m ^2}

Pour l’estimation d’une moyenne:

n = \frac{\epsilon ^2 \cdot \sigma ^2}{m ^2}

Pour la comparaison de 2 proportions ou de deux pourcentages P_1, P_2:

L’hypothèse H_0 testée est P_1 = P_2. La formule devient alors:

n = \frac{\Phi ^2}{2 (Arcsin( \sqrt{P_1}) - Arcsin (\sqrt{P_2})) ^2}

Pour la comparaison de deux moyennes \mu_1, \mu_2:

L’hypothèse H_0 testée est \mu_1 = \mu_2. La formule devient alors:

n = \frac{2 \cdot \sigma ^2 \cdot \Phi ^2}{(\mu_1 - \mu_2) ^2}

Avec \Phi = R(\alpha) + R(2 \beta) dans le cadre bilatéral et \Phi = R(2 \alpha) + R(2 \beta) pour celui unilatéral. Nous rappelons que R(\alpha) représente le risque de première espèce soit celui de rejetter à tort H_0 et R(\beta) celui de seconde espèce soit le risque de ne pas détecter une différence significative pourtant bien présente (ce qui est appelé « manque de puissance » dans le jargon du statisticien).

A noter que l’ensemble de ces formules sont fortement dépendantes de la connaissance à priori du résultat que nous souhaitons mesurer en réalité. Certes ce constat est assez paradoxale mais terriblement logique. Il convient alors de déterminer les différentes estimations de P, P_1, P_2, \sigma ^2 au travers d’un travail bibliographique assidu ou bien d’une pré-enquête.

Enfin, un ajustement lié à l’effet du plan d’échantillonnage est nécessaire dans le cas particulier du sondage à plusieurs degrés (et du sondage en grappes), il convient alors de multiplier le n estimé par le nombre de degrés.

\bullet Les différentes méthodes de sondages

1) Sondage aléatoire simple

Il s’agit de la méthode la plus connue. Elle consiste simplement à tirer aléatoirement et sans remise n <<< N unités de sondage (ou selon sa variante, via un pas ou tirage systémique) depuis la population d’étude de taille N.

Pour qu’un sondage aléatoire fournissent des résultats précis il faut agir:

♦ sur la taille n pour qu’il soit assez grand, ce qui représente la condition la plus importante de ce type de sondage

♦ sur le taux de sondage f = \frac{n}{N} pour qu’il soit le plus possible voisin de 1 (à noter que le cas f = 1 équivaut à un recensement)

♦ sur une variance S ^2 la plus faible possible afin de garantir une précision maximale des estimateurs

Le tirage aléatoire simple constitue la base des plans de sondage, nous verrons que les familles présentées ci-dessous consiste à « craqueler » l’information pour finalement appliquer cette méthode sur les sous-groupes constitués.

Notons également qu’une autre version du tirage aléatoire simple existe, celle avec remise permettant notamment de garantir que sur de grands échantillons (loi des grands nombres) la distribution de l’estimateur suit une loi normale, permettant ainsi de calculer plus facilement et robustement les intervalles de confiance.

2) Sondage stratifié

Les données peuvent être hétérogènes et présenter des particularités propres à des groupes formés en son sein. Le sondage stratifié permet d’adapter le plan de sondage à ces spécificités en composant H groupes homogènes de la population d’étude et en tirant au sein de ces différents groupes l’échantillon afin d’améliorer la précision des estimateurs, les rendant moins dépendant du hasard.

La formule de décomposition de la variance, qui est une somme entre la variance intra-groupe et la variance inter-groupe, prend une importance particulière pour ce type de sondage dont l’objectif est donc de consituer des groupes tels que la dispersion inter-groupe soit la plus grande possible et la dispersion intra-groupe la plus faible possible.

Les strates se font à partir des variables auxiliaires adaptées au type d’étude souhaitée. Le coût n’est pas forcément plus important que celui d’un sondage aléatoire simple, à condition de disposer d’une base de données disposant de l’information auxiliaire à utiliser. Ce coût reste même négligeable par rapport aux autres coûts.

♦ Le sondage est dit à allocation proportionnelle si,

\forall h \in H, \frac{n_h}{n} = \frac{N_h}{N}

, qui offre une meilleure précision que le sondage aléatoire simple.

♦ Le sondage est dit auto-pondéré si,

\forall h \in H, \frac{n_h}{N_h} = \frac{n}{N} = Cste = f

, qui offre des pondérations identiques pour chaque observation tirée.

3) Sondage à plusieurs degrés

Le sondage aléatoire simple souffre de deux prinpales contraintes: le coût de déplacement de l’enquêteur et le besoin de disposer d’une base complète afin d’assurer de la représentativité du tirage. Une solution envisageable reste le sondage à plusieurs degrés.

Il consiste à construire des groupes d’individus disjoints et dont la réunion est égal à la population totale. Il suffit alors de tirer par sondage aléatoire simple (tirage aléatoire ou tirage systématique) un certain nombre de groupes d’individus qui devient alors une unité d’échantillonnage appelée unité primaire (UP) puis tirer aléatoirement dans chacun de ces groupes selectionnés, ce qui constitue les unités secondaires (nous décrivons ici un sondage à deux degrés, si nous continuons le processus nous rajoutons des degrés au fur et à mesure). A noter qu’une bonne construction d’unités primaires a une logique tout à fait contraire à celle qui prévalait lors du découpage d’une population en strates pour le sondage stratifié. Ensuite, plus nous ajoutons de degrés et plus la précision s’en trouve diminuée.

Le gain, qui permet de pallier aux problèmes du sondage aléatoire simple, réside dans le fait que le sondage à plusieurs degrés nécessite uniquement l’information sur les unités primaires tirées. Le sondage à plusieurs degrés restent moins précis que le sondage aléatoire simple car la constitution des unités primaires implique une homogénéité biaisée dû à l’exclusion de certaines (et donc de groupes d’individus complets) lors du tirage. Le tirage nécessite l’information auxiliaire (appartenance à telle ou telle unité primaire) et représente le seul sondage dont elle apporte une amélioration des estimateurs.

L’enjeu du sondage à plusieurs degrés est de consituer des unités primaires de tailles voisines, aussi faibles que possible et de comportement « moyen » semblable, ce qui revient à regrouper au sein d’une même unité primaire des individus très différents car, selon la formule de décomposition de la variance, la dispersion est concentrée sur la partie « intre-UP » et non la partie « inter-UP ».

Le plan de sondage est dit auto-pondéré si les taux de sondage au second degré sont constants, c’est à dire si nous choisissons des unités secondaires de taille proportionnel à la taille des unités primaires correspondantes, impliquant que tous les individus ont même poids de sondage.

Un cas particulier de sondage à plusieurs degrés bien connu est celui du sondage en grappes qui consiste à enquêter de manière exhaustive l’avant dernière unité d’échantillonnage tirée substituant le processus de tirage aléatoire, devenant ainsi ce qui est appelé « des grappes ». Ce type de sondage subit une perte de précision du fait de la similarité entre individu d’une même unité primaire, chaque tirage amène ainsi son effet de grappe car les unités définies aux degrés ultérieurs ont une certaine ressemblance vis à vis de la variable d’intérêt.

L’objectif pour un bon sondage en grappes devient alors:

♦ des grappes hétérogènes permettant de capter au maximum les différentes particularités de la population d’étude,

♦ des grappes de tailles faibles,

♦ des grappes de tailles voisines, nous retrouvons là l’une des particularités communes aux sondages à plusieurs degrés,

♦ un maximum de grappes à tirer, une autre particularités des sondages à plusieurs degrés.

Le dernier point ne dépend pas vraiment du sondeur, mais plus du budget qui lui est alloué.

4) Sondage à probabilités inégales

Le sondage à probabilités inégales concerne le cas où un individu a plus de change d’être tiré qu’un autre. Il existe alors un très grand nombre de méthodes de tirage visant à assurer la représentativité de l’échantillon tiré à partir de la population d’étude. L’objectif de ces outils étant de construire un plan de sondage pour lequel le tirage de l’echantillon offrira à posteriori à chaque individu une probabilité de sélection égale.

Parmi les algorithmes de tirage, les plus connus nous trouvons:

♦ Le tirage systématique sur un fichier de probabilités cummulées. Les étapes sont:

– considérer l’ensemble des individus i \in [1, N] et leur pobabilité d’inclusion respecive P_i \in ]0, 1],

– cumuler les probabilités d’inclusion (latex]P_i[/latex] itérativement de la manière suivant: P_1, P_1 + P_2, \cdots, \sum_{i = 1} ^N P_i,

– tirer un nombre aléatoire U entre 0 et 1,

– retenir les plages d’individus incluses entre deux cumuls de probabilités d’inclusions adjacents tel que (avec A_0 = 0):

A_{k - 1} = \sum_{i = 1} ^{k - 1} P_i \leq U + (k - 1) < \sum_{i = 1} ^k P_i = A_k

♦ La méthode de Rao-Hartley-Cochran (RHC), qui reste une l’un des outils les plus simple et efficace. Les étapes sont:

– trier les individus dans un ordre aléatoire,

– découper en n groupes successifs et de taille identique \frac{n}{N} (le découpage peut imposer que certains groupes aient un peu plus d’unité de sondage mais ça n’alterne pas la qualité de la l’algorithme),

– tirer proprotionnellement à la taille les unités de sondage dans chaque groupe, en générant une variable aléatoire U \in [0,1] uniforme,

– retenir l’échantillon d’individus dont la variable U est comprise entre les deux cumuls d’effectifs (noté T) de deux groupes adjacents tel que,

\sum_{j = 0} ^{i - 1} \frac{T_j}{T} < U \leq \sum_{j = 1} ^i \frac{T_j}{T}

La méthode RHC est l’une des meilleures si les effectifs et la variable d’intérêt sont à peu prés proportionnels, de plus le calcul de l’estimateur de la variance sans biais en devient simple. A noter que cet outil n’est pas à proprement parlé une technique d’échantillonnage globalement proportionnelle à la taille.

♦ L’algorithme de Sunter, qui reste très simple et qui permet le calcul des probabilités d’inclusion doubles et par conséquent un calcul d’estimateur des variances sans biais, les étapes sont:

– trier les individus selon la valeur décroissante de leur taille et attribuer un rang,

– tirer un individu selon un nombre aléatoire U_i \in [0, 1], s’il est inférieur à la nouvelle probabilité d’inclusion alors il est conservé, il faut réitérer cette étape jusqu’à ce que le taille de groupe soit égale à la taille de l’échantillon.

Cet algorithme permet toujours un tirage de taille fixe ainsi que le calcul des probabilités d’inclusion de tout individu et des probabilités d’inclusion doubles (qui en plus conservent certaines propriétés facilitant le calcul de l’estimateur de la variance sans biais).

♦ Le tirage de Poisson ou méthode de Bernouilli si les probabilités d’inclusion sont constantes, cet outil reste peu efficace mais simple à mettre en palce. Les étapes sont:

– attribuer aux individus une probabilité d’inclusion,

– tirer un nombre aléatoire U_i indépendant des individus,

– si ce nombre est inférieur à la probabilité d’inclusion fixée P_i alors il faut le retenir.

Cet outil construit des échantillons de taille variable.

Dans le cadre du tirage à probabilités inégales, il est fréquent de vouloir distinguer plusieurs catégories d’individus, le but étant alors de pouvoir leur attribuer une probabilité d’inclusion en fonction de ces catégories. Le tirage en deux phases avec post-stratification permet d’y procéder. Les étapes sont alors:

– tirer selon un plan complexe (sondage stratifié ou encore à plusieurs degrés),

– ce qui permet de construire des sous-groupes,

– tirage aléatoire au sein des différents sous-groupes.

La méthode permet d’assurer un taux de sondage constant. Elle reste utile quand on ne dispose pas de l’information auxiliaire ou encore que l’étude est soumise à des coûts de collecte excessifs pour y avoir accés. Ainsi la première phase désigne les individus sur lesquels il faut reccueillir l’information auxiliaire (et peut également servir à contrôler la précision des estimateurs obtenus) et permet de tirer l’échantillon lors de la seconde phase.

Une autre situation peut être rencontrée, celle où l’unité d’observation diffère de l’unité d’échantillonnage (mais peut encore être enquêtée), cela implique que les probabilités d’inclusion sont incalculables pour certaines observations. Le sondage indirect est alors utilisé au travers de la méthode partage des poids qui permet une pondération correct mais pas forcément de façon optimal. La méthode est surtout utilisée dans le cas de bases de sondage multiple, c’est à dire qu’il faut passer par plusieurs bases afin de pouvoir reconstituer la population d’étude et les informations nécessaires.

Enfin, le tirage à deux phases généralisé qui est une méthode spécifique de tirage à probabilités inégales. L’algorithme est:

– tirer un échantillon de la population entière selon un pan quelconque,

– construire un système de probabilités de sélection dans l’échantillon et tirer un nombre d’individus dans cet échantillon selon les probabilités d’inclusion de 1.

Il s’agit en réalité de deux tirages successifs emboîtés appelés première et seconde phases. La méthode permet notamment de résoudre les problèmes pratiques de détermination de la taille de l’échantillon dans des tirages à plusieurs degrés.

5) Sondage équilibré

Le sondage équilibré (ou échantillonnage équilibré) porte également le nom de sondage (ou échantillonnage) contrôlé. L’objectif est d’estimer de façon exact un paramètre connu à priori sur l’ensemble de la population et formée à partir d’une variable auxiliaire particulièrement corrélée à la variable d’intérêt, gagnant ainsi en représentativité de l’échantillon sélectionné.

Deux principaux algorithmes se dégagent de ce type de sondage:

♦ L’algorithme à probabilités égales, qui reste le plus simple notamment si nous sommes en présence d’une unique variable auxiliaire.

Il faut ordonner la variable à estimer par ordre croissant et construire deux sous-population en fonction de l’estimateur souhaité. Si la valeur de l’individu est supérieure ou égale à l’estimateur alors il est classé dans un groupe, sinon dans l’autre. Pour l’individu suivant, il faut réestimer le paramètre et appliquer le même processus. La méthode reste réalisable si l’échantillon est suffisamment grand et si la configuration des données présente une certaine régularité.

♦ L’algorithme CUBE, qui fait partie des méthodes rejectives.

Il consiste à tirer un échantillon selon une méthode quelconque, calculer l’estimateur d’intérêt et voir s’il est inférieur à la valeur réelle à seuil fixé. Si c’est le cas alors l’échantillon est sélectionné, sinon il faut en retirer un autre et procéder au même contrôle jusqu’à trouver l’échantillon satisfaisant .L’échantillon retenu est ainsi considéré comme équilibré. Il est évident que cet algorithme est à la fois le plus performant mais également le plus coûteux, il s’agit en réalité d’explorer petit à petit l’ensemble des combinaisons d’individus satisfaisant les contraintes de représentativité recherchées.  De plus, il met en évidence le non contrôle des probabilités de sélection. L’échantillon ainsi tiré s’apparente à la réalisation de quotas probabilistes, plus performant que sa version empirique car il permet de déterminer une méthode aléatoire et sans biais qui assure la réalisation de certains quotas marginaux, ce qui représente un gain de qualité.

6) Sondage empirique

Le sondage empirique s’oppose assez logiquement au sondage probabiliste car il ne permet plus de déterminer à priori quelle est la probabilité qu’a chaque individu de la population d’étude d’appartenir à l’échantillon. Il est souvent utilisé en l’absence de base de sondage. L’objectif est alors de se rapprocher de le plus possible d’un tirage aléatoire. Les principaux attraits du sondage empirique sont sa rapiditié de mise en oeuvre et son coût d’importance moindre. Plusieurs méthodes existent:

♦ La méthode des quotas, qui est la plus fréquente, où il s’agit de faire en sorte que la structure du sous-échantillon soit exactement la même que celle de la population d’étude. Elle se démarque également par le fait que la précision des estimateurs n’est alors plus calculable de part son statut de méthode empirique et non probabiliste. Le procédé est assez simple, il s’agit de laisser l’enquêteur sélectionner lui-même les unités statistiques sous contrainte des différents quotas fixés au préalable. En général, la méthode des quotas est réservée aux petits échantillons quand le sondage probabiliste est adapté aux échantillons de taille plus importante.

La méthode n’est acceptable que si les principales variables expliquant la variable d’intérêt sont pris en compte comme critère de quotas sans risque d’occasionner un biais de sélection.

Deux types de quotas peuvent être utilisés, les quotas marginaux si l’information croisé entre variables auxiliaires n’est pas disponible, quotas croisés si c’est le cas.

♦ La méthode des itinéraires, semblable à celle des quotas, diffère par le fait que l’enquêteur se voit fixé une contrainte de plus au dépend d’un coût de préparation supplémentaire. Ainsi, il se voit dans l’obligation de prendre en compte le lieu géographique afin de pouvoir sélectionner une unité statistique.

♦ La méthode des unités-types, la plus empirique des méthodes de cette famille de sondage. Elle consiste à sélectionner un individu « moyen » représentatif d’une catégorie complète d’individu (nommé unité-type).

♦ L’échantillonnage volontaire, peu onéreux mais ne permettant pas d’avoir de contrôle sur la population qui ne se portera pas volontaire à l’enquête.

\bullet Application informatique:

Procédure SAS: http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#surveyselect_toc.htm

Package et fonction R: https://cran.r-project.org/web/packages/survey/survey.pdf

\bullet Bibliographie:

– Les techniques de Sondage de Pascal Ardilly

– Sampling with unequal probabilites de K. R. W. Brewer et M. Hanif

– Probabilité, analyse des données et Statistiques de Gilbert Saporta

– Statistique, dictionnaire encyclopédique de Yadolah Dodge

– Médecine tropicale de Marc Gentilini, Eric Caumes, Martin Danis, Dominique Richard-Lenoble, Pierre Bégué, Jean-Etienne Touzé et Dominique Kerouédan

 – Introduction au nombre de sujets nécessaires de Yohann Foucher

– Calcul du nombre de sujets nécessaires (NSN) pour un essai de supériorité d’Olivier Chassany