Comparaison de certains estimateurs de variance positifs pour le modèle d’estimation sur petits domaines Fay-Herriot 6. Résultats de la simulation et analyseComparaison de certains estimateurs de variance positifs pour le modèle d’estimation sur petits domaines Fay-Herriot 6. Résultats de la simulation et analyse

• Précédent
• Suivant

6.1 Distribution Monte Carlo des estimateurs de variance

Le tableau 6.1 montre que l’estimateur de variance REML présente le biais le plus faible $\left({\sigma }_v^2=1\right)$ et la variance la plus élevée. L’efficacité plus faible du REML pourrait être attribuable au fait qu’il ne s’agit pas d’une fonction lisse des données causée par sa définition fractionnée (3.1). L’estimateur MIX hérite d’une partie de cette efficacité faible. Les autres estimateurs de variance ont une variabilité plus faible et un biais positif plus élevé, mais l’espérance conditionnelle des estimateurs AM.YL et AR.YL étant donné que ${\widehat{\sigma }}_{\text{vREML}}^{\text{2}}=0$ est proche de zéro. Le biais inconditionnel de l’estimateur AM.LL est plus élevé que celui du MIX. Selon la définition de l’estimateur MIX, les biais conditionnels des estimateurs MIX et AM.LL coïncident. En outre, le MIX converge plus rapidement que les autres estimateurs. Par exemple, étant donné la distribution des probabilités sur les 10 000 estimations de la variance où $m=45,$ nous avons calculé la probabilité que les estimations se situent dans un intervalle contenant ${\sigma }_v^2=1.$ La probabilité que les estimations se situent entre 0,6 et 1,4 est de 0,47 pour le MIX et de 0,16 pour l’AM.YL. Par ailleurs, la probabilité que les estimations soient inférieures à 0,2 est de 0,05 pour le MIX et de 0,53 pour l’AM.YL.

6.2 EQM réelle de l’EBLUP, biais relatif moyen et racine de l’EQM relative moyenne des estimateurs de l’EQM

Tous les estimateurs de variance sont convergents et asymptotiquement normaux, la variance convergeant au même rythme. Leurs biais sont différents : ceux des estimateurs REML, AR.YL et MIX sont de l’ordre  $o\left(1/m\right),$ tandis que ceux des estimateurs AM.LL et AM.YL sont de l’ordre  $O\left(1/m\right).$ Le biais inhérent aux trois dernières méthodes ont un impact sur l’estimation de l’EQM de l’EBLUP, même pour un nombre modéré de domaines.

Pour $m=100,$ les tableaux 6.2a et 6.2b montrent que l’EQM de l’EBLUP diminue à mesure que ${\sigma }_v^2/{\psi }_i$ augmente, et que cette relation se maintient quel que soit le nombre de domaines. Nous observons que l’EQM de ${\widehat{\theta }}_i$ sous les estimateurs de variance REML et MIX est légèrement plus élevée que le reste des EQM en raison de la plus grande variabilité inhérente à ces estimateurs de variance. Le tableau 6.2a présente les résultats pour l’estimateur de l’EQM par linéarisation de Taylor et les deux estimateurs paramétriques de l’EQM sous les estimateurs de variance REML, AM.LL, AR.YL et AM.YL. Le tableau 6.2b présente les résultats pour les estimateurs suivants de l’EQM sous l’estimation de variance MIX : RB_Y1 défini en (4.3), RB_Y2 défini en (4.2), M_et_coll défini en (4.5), BP EQM et BP EQM naïf. Parmi les estimateurs de l’EQM de Taylor, RB_Y1 et M_et_coll sous la méthode MIX présentent le biais le plus faible. Parmi les estimateurs bootstrap de l’EQM, BP sous MIX et BP naïf sous AR.YL présentent le biais le plus faible. Quant à la racine de l’erreur quadratique moyenne relative (REQMR) des estimateurs de l’EQM, elle diminue à mesure que ${\sigma }_v^2/{\psi }_i$ augmente. Les différences entre l’estimateur RB_Y2 de l’EQM sous le MIX et l’estimateur de l’EQM de Taylor sous l’AM.YL semblent faibles mais consistantes. Alors que le RB_Y1, le M_et_coll et les estimateurs naïfs de l’EQM sous la méthode MIX présentent un BRM plus faible que le RB_Y2 sous la même méthode, et que le BRM de l’estimateur de Taylor et de l’estimateur BP naïf sous la méthode AR.YL est plus faible que celui du RB_Y2 sous la méthode MIX, c’est le contraire pour la REQMR. Cela s’explique en partie par le biais conditionnel négatif extrême de ces estimateurs de l’EQM (c’est-à-dire les estimateurs RB_Y1 et M_et_coll sous la méthode MIX et les estimateurs de Taylor et BP naïf sous la méthode AR.YL), comme le montre le tableau 6.3. Même pour $m=100,$ une proportion relativement élevée (16 %) des populations donnent ${\widehat{\sigma }}_{v\text{REML}}^2=0$ et, dans ces populations, les estimations obtenues au moyen de la plupart des méthodes d’estimation de variance et la plupart des estimateurs de l’EQM sont les plus inférieures à la valeur réelle. C’est-à-dire que, pour ces estimateurs de l’EQM, les estimateurs conditionnels ne donnent pas de bons résultats. L’estimateur BP EQM semble corriger pour le biais de façon satisfaisante, mais il est plus variable que le BP EQM naïf. Lorsque nous incluons le BRM, la REQMR et le ${\text{BRM}}_{ℂ}$ dans l’évaluation, c’est le RB_Y2 sous la méthode MIX, suivi de près par l’estimateur BP naïf, qui l’emporte sous la méthode MIX qui semble donner les meilleurs résultats. Nous pourrions donc conclure à la supériorité des estimateurs RB_Y2 et naïf sous la méthode MIX pour $m=100,$ ce qui est un nombre modéré de domaines pour les données de ce genre.

Les tableaux 6.4a et 6.4b ci-dessous présentent les résultats pour $m=45$ avec 9 domaines par ${\sigma }_v^2/{\psi }_i.$ L’estimateur AM.YL donne des EQM plus petites que le MIX, les différences ne dépassant pas 2 %. Le biais des estimateurs de variance augmente à mesure que le nombre de domaines diminue, ce qui a un impact sur les estimateurs de l’EQM. En effet, les BRM de tous les estimateurs de l’EQM ont augmenté. En particulier, le BRM des estimateurs de l’EQM de Taylor sous l’estimation de variance YL et LL et le BRM de l’estimateur RB_Y2 ont augmenté de 100 % par rapport au BRM avec 100 domaines. En ce qui concerne la REQMR, l’estimateur de l’EQM de Taylor sous AM.YL a une REQMR légèrement inférieure à celle du RB_Y2 sous la méthode MIX pour une ${\sigma }_v^2/{\psi }_i$ très faible. En général, la variabilité (en termes de REQMR) du RB_Y2 est plus faible que celle de l’estimateur de Taylor sous LL et YL et des estimateurs RB_Y1 et M_et_coll. Cela pourrait être attribuable en partie à la sous-estimation des EQM pour les populations avec des estimations REML nulles, dont le pourcentage tourne autour de 30 % lorsque $m=45.$ Le tableau 6.5 est plus instructif à cet égard : étant donné ${\widehat{\sigma }}_{vR\text{EML}}^2=0,$ RB_Y1 et M_et_coll conduisent à une sous-estimation importante.

Compte tenu du BRM, de la REQMR et du ${\text{BRM}}_{ℂ}$ des estimateurs de l’EQM, c’est l’estimateur BP naïf de l’EQM sous la méthode MIX qui donne les meilleurs résultats pour les ${\sigma }_v^2/{\psi }_i$ plus importantes. Le tableau 6.6 présente la moyenne des mesures de performance, calculée sur les cinq groupes de variances d’échantillonnage, pour les trois estimateurs de l’EQM de Taylor sous la méthode MIX avec des données du modèle décrit en 5.1, mais avec trois valeurs différentes de ${\sigma }_v^2.$ L’estimateur RB_Y2 donne de meilleurs résultats quand ${\sigma }_v^2=1,$ mais lorsque ${\sigma }_v^2$ diminue, c’est l’estimateur de l’EQM de M_et_coll qui prend le dessus, précisément parce qu’il a été construit en partant du principe que ${\sigma }_v^2$ est d’environ zéro.

Les tableaux 6.7a et 6.7b ci-dessous montrent les résultats pour $m=15$ domaines avec 3 domaines par ${\sigma }_v^2/{\psi }_i.$ Les EQM ne varient pas de plus de 5 %, quelle que soit la méthode d’estimation de variance utilisée.

Il n’existe pas de relation monotone entre le BRM ou la REQMR et ${\sigma }_v^2/{\psi }_i,$ ce qui pourrait indiquer que l’approximation d’ordre deux pour estimer l’EQM est médiocre, quelle que soit la méthode d’estimation de variance utilisée. Les BRM des estimateurs de l’EQM de Taylor sous les méthodes d’estimation de variance de LL et YL sont trop élevés, et il en va de même pour la REQMR. L’estimateur RB_Y2 sous la méthode MIX ne s’en sort pas très bien non plus. La raison de ce résultat est claire : le pourcentage élevé d’estimations REML nulles (43 %) suggère que le MIX coïncide avec l’AM.LL pour les populations REML nulles. Ainsi, le MIX a un biais positif pour $m=15,$ et le RB_Y2 ne tient pas compte de ce biais. Le RB_Y1 tient compte du biais dans le MIX, mais l’estimateur du biais n’est pas très précis pour $m=15.$ L’estimateur de l’EQM de M_et_coll coïncide presque avec le BRM et la REQMR de l’estimateur de l’EQM de Taylor sous l’estimation de variance REML, car ils sont égaux par définition lorsque ${\widehat{\sigma }}_{v\text{REML}}^2=0.$ Le ${\text{BRM}}_{ℂ}$ des trois estimateurs de l’EQM de Taylor sous la méthode MIX laisse à désirer. Compte tenu de toutes les mesures de performance, les estimateurs bootstrap de l’EQM donnent de meilleurs résultats que les estimateurs de l’EQM de Taylor. Pour $m=15$ domaines avec 3 domaines par ${\sigma }_v^2/{\psi }_i,$ l’estimateur BP sous la méthode MIX donne les meilleurs résultats, l’estimateur BP naïf sous AR.YL et AM.YL venant en deuxième place.

En résumé, sous le modèle de Fay-Herriot avec une ${\sigma }_v^2$ positive, MIX et AR.YL sont les seuls estimateurs de variance positifs à l’étude qui présentent un biais asymptotique négligeable. Le biais asymptotique des estimateurs de variance AM.YL et LL est plus grand. En revanche, notre simulation a démontré que, pour un nombre modéré de domaines et pour les populations qui donnent des estimations REML nulles, les deux estimateurs de variance de YL présentaient un biais négatif et produisaient des EBLUP proches de l’estimateur synthétique de la moyenne. Par contre, le MIX, qui combine les estimateurs AM.LL et REML, ne présentait qu’un biais légèrement négatif dans ces populations. De plus, la distribution inconditionnelle du MIX approchait de la normalité beaucoup plus rapidement que celle des autres estimateurs de variance.

$${}_{}$$ En ce qui concerne l’estimateur de l’EQM de l’EBLUP, il était beaucoup plus précis que l’estimateur direct, sous toutes les méthodes d’estimation de variance examinées ici, même pour un petit nombre de domaines. Les estimateurs de variance AM.LL, AM.YL et AR.YL étaient tous moins variables que le REML et le MIX. L’impact sur l’EQM de l’EBLUP était minime, car il y avait peu de différences entre les EQM pour le même rapport signal/bruit. Ces différences s’accentuaient à mesure que le nombre de domaines ou le rapport signal/bruit diminuait. Ainsi, pour un rapport signal/bruit extrêmement faible, l’EQM sous la méthode MIX pourrait être un peu plus grande que sous l’estimateur de variance AM.YL.

Sous la méthode MIX d’estimation de variance, nous avons comparé trois estimateurs de l’EQM de type Taylor et deux estimateurs de l’EQM bootstrap. Les trois estimateurs de l’EQM de Taylor sous la méthode MIX (RB_Y1, RB_Y2 et M_et_coll) sont sans biais jusqu’à l’ordre deux. Les estimateurs de l’EQM de type Taylor sous LL et YL sont eux aussi sans biais jusqu’à l’ordre deux. Les estimateurs RB_Y1, AM.LL et AM.YL peuvent produire des estimations négatives de l’EQM.

L’estimateur de l’EQM de Taylor sous la méthode REML d’estimation de variance et l’estimateur M_et_coll sous la méthode MIX coïncidant par définition, les différences entre leurs mesures de performance sont négligeables (leurs EQM réelles sont différentes; dans notre étude cependant, pour $m=100,$ l’estimateur MIX coïncidait avec le REML 84 % du temps). Pour un nombre modéré de domaines, qui pourrait être $m=45\text{ ou }100$ pour ces données, et pour les populations qui donnent des estimations REML nulles, les deux estimateurs de l’EQM de Taylor sous le REML et les estimateurs de l’EQM de M_et_coll ne tiennent pas compte de la variation attribuable à l’estimation de ${\sigma }_v^2,$ ce qui se reflète dans leur ${\text{BRM}}_{ℂ},$ très négatif, qui est inférieur à -60 % pour les rapports signal/bruit plus petits. Par ailleurs, l’estimateur RB_Y1 tient compte de la variation due à l’estimation de ${\sigma }_v^2,$ mais son ${\text{BRM}}_{ℂ}$ est lui aussi très négatif. En effet, le RB_Y1 est un estimateur de l’EQM fractionné qui, pour les populations où ${\widehat{\sigma }}_{v\text{REML}}^2=0,$ soustrait un facteur du biais inconditionnel de l’AM.LL, qui est toujours positif. Une meilleure formule pour un estimateur de l’EQM fractionné serait d’utiliser un estimateur du biais conditionnel $E\left({\widehat{\sigma }}_v^2/{\widehat{\sigma }}_{v\text{REML}}^2=0\right).$ En fait, même pour un nombre modéré de domaines $\left(m=100\right),$ le tableau 6.1 montre que l’estimateur MIX a un biais inconditionnel de 49 %, mais un biais conditionnel de -37 %.

L’estimateur BP de l’EQM sous les méthodes AR.YL et MIX était bien corrigé pour le biais, mais la variance en souffrait. Le bootstrap naïf semble être l’estimateur de l’EQM qui donne les meilleurs résultats, et ces résultats sont encore meilleurs sous l’estimation de variance MIX, même lorsque les trois mesures (BRM, ${\text{BRM}}_{ℂ}$ et REQMR) sont prises en compte. Nous avons constaté que, pour un nombre modéré de domaines, le RB_Y2 était l’estimateur de Taylor sous la méthode MIX qui présentait la REQMR la plus faible. Par ailleurs, l’estimateur de M_et_coll est le plus fiable lorsque la variance sous-jacente réelle ${\sigma }_v^2$ est très faible : dans ce cas-ci, M_et_coll est effectivement l’estimateur de l’EQM de l’estimateur synthétique de la moyenne de petit domaine. Nous ne recommandons pas de compter sur l’approximation d’ordre deux de l’EQM lorsque $m$ est petit : l’approximation (2.6) de l’EQM ne tient pas nécessairement, les mesures de performance obtenues dans notre étude sont très instables, et elles peuvent varier d’un ensemble de données à l’autre.

En conclusion, sous l’hypothèse ${\sigma }_v^2>0,$ les performances relatives des estimateurs de variance positifs comparés dépendent de la taille de ${\sigma }_v^2,$ du rapport signal/bruit, du nombre de domaines et de la fonction objectif. Pour un nombre modéré de domaines, l’estimateur de variance MIX semblait donner de meilleurs résultats que les estimateurs LL et YL dans cette étude. Sous la méthode MIX, l’estimateur BP naïf de l’EQM présentait la combinaison ${\text{BRM}}_{ℂ}$ et REQMR la plus faible. L’estimateur de l’EQM de M_et_coll sous l’estimateur de variance MIX donnait des résultats légèrement meilleurs que le RB_Y1 lorsque la variance ${\sigma }_v^2$ sous-jacente était très petite. Cependant, le pourcentage de REML nuls produits sous le modèle de simulation montre qu’un résultat de ${\widehat{\sigma }}_{v\text{REML}}^2=0$ et/ou que des tests d’hypothèses négatifs ne signifient pas nécessairement que la variance ${\sigma }_v^2$ est assez faible pour que l’on puisse se fier à l’estimateur de M_et_coll. En l’absence d’autres renseignements, l’estimateur BP naïf sous la méthode MIX semble donner de meilleurs résultats.

Remerciements

Les auteurs désirent remercier le professeur J.N.K. Rao de l’Université Carleton pour ses commentaires utiles ainsi que Victor Estevao de Statistique Canada, qui a développé un algorithme de maximisation de grille spécialement pour ce projet. Ils tiennent également à remercier les membres du comité de revue pour leur évaluation consciencieuse de cet article et pour leurs suggestions d’amélioration.

Annexe A

Preuve du théorème 4.1

La variance asymptotique de ${\widehat{\sigma }}_{v\text{MIX}}^2$ est donnée par : $\overline{V}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)=\underset{m\to \infty }{\lim }E{\left({\widehat{\sigma }}_{v\text{MIX}}^2-{\sigma }_v^2\right)}^2$

Nous montrons que $E{\left({\widehat{\sigma }}_{v\text{MIX}}^2-{\sigma }_v^2\right)}^2\le E{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^2+o\left(1/m\right)$ à mesure que $m\to \infty .$

$$\begin{array}{ll}E{\left({\widehat{\sigma }}_{v\text{MIX}}^2-{\sigma }_v^2\right)}^2\hfill & ={\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2>0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^{2}dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2}dP}\hfill \\ \hfill & \begin{array}{l}\le {\displaystyle \underset{\Omega }{\int }{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^{2}dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2}dP}=E{\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)}^2\text{(A}\text{.1)}\\ +o\left(\frac{1}{m}\right).\end{array}\hfill \end{array}$$

En effet, par les inégalités de Holder et Minkowski, avec $1<p<\infty ,1/p+1/q=1,$ et si $X\equiv {\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^2=O_p\left(1/m\right)$ et l’indicateur $I\left({\widehat{\sigma }}_{\text{vREML}}^{\text{2}}=0\right)$ des populations avec ${\widehat{\sigma }}_{v\text{REML}}^2=0,$ nous avons :

$$\begin{array}{ll}{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2<0\right\}}{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2}dP}\hfill & \le {\left({\displaystyle \underset{\Omega }{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^{2p}dP}\right)}^{1/p}.{\left(P\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}\right)}^{1/q}\hfill \\ \hfill & ={\left(O\left(\frac{1}{m^p}\right)\right)}^{1/p}.{\left(o\left(1\right)\right)}^{1/q}=o\left(\frac{1}{m}\right),\hfill \end{array}(\text{A}.2)$$

puisque ${\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)}^2$ est uniformément borné et ${\widehat{\sigma }}_{v\text{REML}}^2\stackrel{P}{{\displaystyle \to }}{\sigma }_v^2>0.$ Il est à noter que les estimateurs AM.LL et REML de ${\sigma }_v^2$ sont uniformément bornés en conséquence de leur convergence presque certaine vers ${\sigma }_v^2$ (voir, par exemple, Yuan et Jennrich 1998).

Preuve du théorème 4.2

Nous désignons par ${\widehat{\sigma }}_{v\text{ML}}^2$ l’estimateur de variance par le maximum de vraisemblance.

Nous montrons d’abord que ${\widehat{\sigma }}_{v\text{REML}}^2-{\widehat{\sigma }}_{v\text{ML}}^2=O_p\left(1/m\right).$ Soit $G_{*}\left({\sigma }_v^2\right)=\partial \log \left(L_{*}\right)/\partial {\sigma }_v^2=0$ l’équation d’estimation qui donne l’estimateur de variance *. L’équation (3.4) implique que :

$$G_{\text{AM}\text{.LL}}\left({\sigma }_v^2\right)-G_{\text{ML}}\left({\sigma }_v^2\right)=\partial \log {\sigma }_v^2/\partial {\sigma }_v^2=\frac{1}{m{\sigma }_v^2}=O\left(\frac{1}{m}\right).(\text{A}.3)$$

Avec $G_{\text{ML}}^{\prime }(\cdot )\triangleq \left(\partial G_{\text{ML}}/\partial {\sigma }_v^2\right)(\cdot )$ et $G_{\text{ML}}^{\prime \text{}\prime }(\cdot )\triangleq \left(\partial G_{\text{ML}}^{\prime }/\partial {\sigma }_v^2\right)(\cdot ),$ l’équation (A.3) implique que :

$$G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)-G_{\text{AM}\text{.LL}}^{\prime }\left({\sigma }_v^2\right)=O\left(\frac{1}{m}\right).(\text{A}.4)$$

Maintenant, en utilisant l’équation (A.4), la $\sqrt{m}-$ convergence des estimateurs ML et AM.LL de ${\sigma }_v^2,$ le développement en série de Taylor à deux termes de $G_{\text{ML}}(\cdot )\text{ et }{G}_{\text{AM}\text{.LL}}(\cdot )$ à ${\sigma }_v^2$ et $G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)=O\left(1\right)$ à mesure que $m\to \infty ,$ le côté gauche de l’équation en (A.3) est égal à:

$$\begin{array}{l}=G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{ML}}^2-{\sigma }_v^2\right)-G_{\text{AM}\text{.LL}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)+O_p\left(\frac{1}{m}\right)\\ =G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{ML}}^2-{\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2\right)+\left(G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)-G_{\text{AM}\text{.LL}}^{\prime }\left({\sigma }_v^2\right)\right)\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)+O_p\left(\frac{1}{m}\right)\\ =G_{\text{ML}}^{\prime }\left({\sigma }_v^2\right)\left({\widehat{\sigma }}_{v\text{ML}}^2-{\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2\right)+O_p\left(\frac{1}{m^{3/2}}\right)+O_p\left(\frac{1}{m}\right).\end{array}$$

La dernière égalité ci-dessus implique que :

$${\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{ML}}^2=O_p\left(\frac{1}{m}\right)\text{ à mesure que }m\to \infty .(\text{A}.5)$$

De même, nous établissons une relation entre $G_{\text{REML}}\left({\sigma }_v^2\right)$ et $G_{\text{ML}}\left({\sigma }_v^2\right):$ étant donné que $\text{tr}\left(V^{-1}Z{\left(Z^{\prime }{V}^{-1}Z\right)}^{-1}{Z}^{\prime }{V}^{-1}\right)=O\left(1\right)$ découle des conditions 1 à 3 de la section 3 et de l’équation (3.1), nous avons :

$$G_{\text{REML}}\left({\sigma }_v^2\right)-G_{\text{ML}}\left({\sigma }_v^2\right)=\frac{1}{m}\text{tr}\left(V^{-1}Z{\left(Z^{\prime }{V}^{-1}Z\right)}^{-1}{Z}^{\prime }{V}^{-1}\right)=O\left(\frac{1}{m}\right)\text{ à mesure que }m\to \infty ,(\text{A}.6)$$

l’équation (A.6) et le même argument que pour l’estimateur AM.LL impliquent que :

$${\widehat{\sigma }}_{v\text{REML}}^2-{\widehat{\sigma }}_{v\text{ML}}^2=O_p\left(\frac{1}{m}\right)\text{ à mesure que }m\to \infty .(\text{A}.7)$$

Ensemble, les équations (A.5) et (A.7) donnent :

$$\left({\widehat{\sigma }}_{v\text{REML}}^2-{\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2\right)=O_P\left(\frac{1}{m}\right).(\text{A}.8)$$

Nous exprimons maintenant le biais de l’estimateur MIX comme suit :

$$B_{\text{MIX}}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)={\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2>0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\sigma }_v^2\right)dP}.$$

Nous ajoutons et soustrayons ${\displaystyle {\int }_{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)dP}$ du côté droit de l’équation ci-dessus pour obtenir :

$$\begin{array}{ll}{B}_{\text{MIX}}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)\hfill & ={\displaystyle \underset{\Omega }{\int }\left({\widehat{\sigma }}_{v\text{REML}}^2-{\sigma }_v^2\right)dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{REML}}^2\right)dP}\hfill \\ \hfill & =\text{Biais}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{REML}}^2\right)dP}.\hfill \end{array}(\text{A}.9)$$

Puisque ${\widehat{\sigma }}_{v\text{AM}\text{.LL}}^{\text{2}}-{\widehat{\sigma }}_{v\text{REML}}^2$ est uniformément borné, nous appliquons l’inégalité de Holder et Minkowski avec $p=q=2$ et l’équation (A.8) au dernier terme en (A.9) pour obtenir :

$$\begin{array}{ll}{B}_{\text{MIX}}\left({\widehat{\sigma }}_{v\text{MIX}}^2\right)\hfill & =\text{Biais}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+{\left({\displaystyle \underset{\Omega }{\int }{\left({\widehat{\sigma }}_{v\text{AM}\text{.LL}}^2-{\widehat{\sigma }}_{v\text{REML}}^2\right)}^{2}dP}\right)}^{1/2}\cdot P{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}^{1/2}\hfill \\ \hfill & =\text{Biais}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+O\left(\frac{1}{m}\right)\cdot o\left(1\right)=\text{Biais}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)+o\left(\frac{1}{m}\right).\hfill \end{array}(\text{A}.10)$$

Preuve de la remarque 4.2 : ${\text{eqm}}_0$ est sans biais jusqu’à l’ordre deux

$$\begin{array}{ll}E\left({\text{eqm}}_0\right)-\text{EQM}\left({\widehat{\theta }}_i\right)\hfill & ={\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2>0\right\}}{\int }\left(g_{1i}+g_{2i}+2g_{3i}\right)\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP}+{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{g}_{2i}}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP-\text{EQM}\hfill \\ \hfill & \begin{array}{l}=\left[{\displaystyle \underset{\Omega }{\int }\left(g_{1i}+g_{2i}+2g_{3i}\right)\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP}-\text{EQM}\right]\\ +{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }{g}_{2i}}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP-{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }\left(g_{1i}+g_{2i}+2g_{3i}\right)\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP}\end{array}\hfill \\ \hfill & =\left[o\left(\frac{1}{m}\right)\right]-{\displaystyle \underset{\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}}{\int }2g_{3i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)dP},\hfill \end{array}(\text{A}.11)$$

puisque $g_{1i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)=g_{1i}\left(0\right)=0$ en $\left\{{\widehat{\sigma }}_{v\text{REML}}^2=0\right\}$ et $g_{2i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)$ s’annulent en (A.11). Mais

$$g_{3i}\left({\widehat{\sigma }}_{v\text{REML}}^2\right)=g_{3i}\left(0\right)=\frac{\overline{V}\left(0\right)}{{\psi }_i}=O_p\left(\frac{1}{m}\right)$$

et est uniformément borné sous les conditions de régularité données dans la section 2. Le dernier terme en (A.11) est donc aussi un $o\left(1/m\right),$ ce qui rend ${\text{eqm}}_0$ sans biais jusqu’à l’ordre deux.

Annexe B

B.1 Comparaison entre les estimateurs REML et AR.YL au moyen de l’algorithme de score

L’algorithme de score donnait parfois des estimations nulles pour la vraisemblance de l’estimateur AR.YL. En effet, pour les ensembles de données simulés sous le modèle donné dans la section 5, où $m=45$ et ${\sigma }_v^2=1,$ les algorithmes de score REML et AR.YL produisaient 28 % et 26 % d’estimations nulles respectivement. On peut voir pourquoi dans les figures B.1 à B.3 : les vraisemblances correspondent à une seule population générée sous le modèle avec ${\sigma }_v^2=1$ pour lequel ${\widehat{\sigma }}_{v\text{REML}}^2=0.$

La figure B.2 montre que la valeur maximale de la vraisemblance AR.YL est très proche de la limite. Il arrive parfois que l’algorithme de score passe à côté du maximum et donne une valeur nulle. La figure B.3 montre que la vraisemblance AM.LL a une valeur maximale qui se différencie mieux de la limite.

B.2 Traitement des zéros dans l’estimateur bootstrap paramétrique

Pour chaque estimation ${\widehat{\sigma }}_v^2={\widehat{\sigma }}_v^2\left(y^{\left(r\right)}\right),r=1,\dots 10K,$ et chaque méthode d’estimation de variance :

1. Générer un grand nombre B d’effets aléatoires de domaine $v_i^{\left(b\right)}\stackrel{\text{i}\text{.i}\text{.d}\text{.}}{\sim }N\left(0,{\widehat{\sigma }}_v^2\right),b=1,\dots ,B,$ et générer, indépendamment de $v_i^{\left(b\right)},$ des erreurs d’échantillonnage $e_i^{\left(b\right)}\stackrel{\text{i}\text{.i}\text{.d}\text{.}}{\sim }N\left(0,{\psi }_i\right),$ $i=1,\dots ,m,b=1,\dots ,B.$ Générer des données bootstrap $y_i^{\left(b\right)}={\theta }_i^{\left(b\right)}+e_i^{\left(b\right)},$ ${\theta }_i^{\left(b\right)}=x_i^{\prime }\widehat{\beta }+v_i^{\left(b\right)},$ $i=1,\dots ,m.$ Si ${\widehat{\sigma }}_{v\text{REML}}^2\left(y^{\left(r\right)}\right)=0,$ générer $\left(y_i^{\left(b\right)},{\theta }_i^{\left(b\right)}\right),$ $b=1,\dots ,B,$ à partir du modèle synthétique (voir aussi Rao et Molina 2015).
2. Adapter le modèle aux données bootstrap et obtenir ${\widehat{\sigma }}_v^{2\left(b\right)};$ pour l’estimateur MIX, calculer ${\widehat{\sigma }}_{v\text{MIX}}^{2\left(b\right)}={\widehat{\sigma }}_{v\text{REML}}^{2\left(b\right)}\text{ si }{\widehat{\sigma }}_v^{2\left(b\right)}$ est positif et ${\widehat{\sigma }}_{\text{vMIX}}^{\text{2}\left(b\right)}={\widehat{\sigma }}_{v\text{AM}}^{2\left(b\right)}$ autrement.
3. Obtenir ${\widehat{\beta }}^{\left(b\right)},$ l’EBLUP correspondant ${\widehat{\theta }}_i^{\left(b\right)},$ les composantes bootstrap $g_{1i}^{\left(b\right)}=g_{1i}\left({\widehat{\sigma }}_v^2\left(y^{\left(b\right)}\right)\right),$ $g_{2i}^{\left(b\right)}=g_{2i}\left({\widehat{\sigma }}_v^2\left(y^{\left(b\right)}\right)\right)$ et ${\overline{g}}_{ji}^{\text{BP}}=B^{-1}{\displaystyle {\sum }_b{g}_{ji}^{\left(b\right)}},j=1,2.$
4. L’estimateur bootstrap naïf de l’EQM est ${\text{eqm}}_{\text{naive}}=B^{-1}{\displaystyle {\sum }_{b=1}^B{\left({\widehat{\theta }}_i^{\left(b\right)}-{\theta }_i^{\left(b\right)}\right)}^2}.$
5. L’estimateur BP de l’EQM (qui est corrigé du biais (Pfeffermann et Glickman 2004)) est : ${\text{eqm}}_{\text{BP}}\left({\widehat{\theta }}_i\right)=g_{1i}\left({\widehat{\sigma }}_v^2\right)+g_{2i}\left({\widehat{\sigma }}_v^2\right)-{\overline{g}}_{1i}^{\text{BP}}-{\overline{g}}_{2i}^{\text{BP}}+{\text{eqm}}_{\text{naive}}.$
6. Pour calculer le ${\text{BRM}}_{ℂ},$ faire la moyenne de $\left({\text{eqm}}_{\text{BP}}^{\left(r\right)}\left({\widehat{\theta }}_i\right)-\text{EQM}\left({\widehat{\theta }}_i\right)\right)/\text{EQM}\left({\widehat{\theta }}_i\right)$ pour les populations où $\left(r\right)/{\widehat{\sigma }}_{v\text{REML}}^2\left(y^{\left(r\right)}\right)=0$ et faire de même pour le ${\text{BRM}}_{ℂ}$ de ${\text{eqm}}_{\text{naive}}.$

Bibliographie

Chen, S., et Lahiri, P. (2008). On mean squared prediction error estimation in small area estimation problems. Communications in Statistics-Theory and Methods, 37, 1792-1798.

Chen, S., et Lahiri, P. (2011). On the estimation of Mean Squared Prediction Error in small area estimation. Calcutta Statistical Association Bulletin, 63, (Special 7^th Triennial Proceedings Volume), Nos. 249-252.

Cressie, N. (1992). Estimation du maximum de vraisemblance avec contrainte (MVC) dans le lissage des taux de sous-dénombrement du recensement selon l’approche empirique de Baye. Techniques d’enquête, 18, 1, 83-103.

Das, K., Jiang, J. et Rao, J.N.K. (2004). Mean squared error of empirical predictor. The Annals of Statistics, 32, 2, 818-840.

Datta, G., et Lahiri, P. (2000). A unified measure of uncertainty of estimated best linear unbiased predictors in small area estimation problems. Statistica Sinica, 10, 613-627.

Estevao, V. (2014). Grid optimization algorithm for maximum likelihood. Rapport interne, Division de la recherche et de l’innovation en statistique (DRIS), Statistique Canada.

Fay, R.E., et Herriot, R.A. (1979). Estimation of income from small places: An application of James-Stein Procedures to census data. Journal of the American Statistical Association, 74, 269-277.

Lahiri, P., et Li, H. (2009). Generalized maximum likelihood method in linear mixed models with an application in small area estimation. Dans Proceedings of the Federal Committee on Statistical Methodology Research Conference, disponible au http//www.fcsm.gov/events/papers2009.html.

Lahiri, P., et Pramanik, S. (2011). Discussion of “Estimating random effects via adjustment for density maximization” par C. Morris et R. Tang. Statistical Science, 26, 2, 291-295.

Li, H., et Lahiri, P. (2011). An adjusted maximum likelihood method for solving small area estimation problems. Journal of Multivariate Analysis, 101, 882-892.

Molina, I., Rao, J.N.K. et Datta, G.S. (2015). Estimation sur petits domaines sous un modèle de Fay-Herriot avec test préliminaire pour la présence d’effets aléatoires de domaine. Techniques d’enquête, 41, 1, 1-20.

Morris, C.N. (2006). Mixed model prediction and small area estimation (with discussions). Test, 15, 72-76.

Pfeffermann, D., et Glickman, H. (2004). Mean squarred error approximation in small area estimation by use of parametric and non-parametric bootstrap. Proceedings of the American Statistical Association, Section on Survey Research Methods, Alexandria, VA. 4167-78.

Rao, J.N.K. (2003). Small Area Estimation. New York : John Wiley & Sons, Inc.

Rao, J.N.K., et Molina, I. (2015). Small Area Estimation, second edition. New York : John Wiley & Sons, Inc.

Rubin-Bleuer, S., et Schiopu-Kratina, I. (2005). On the two-phase framework for joint model and design-based inference. The Annals of Statistics, 33, 6, 2789-2810.

Rubin-Bleuer, S., et You, Y. (2012). A positive variance estimator for the Fay-Herriot small area model. SRID-2012-009E, Division de la recherche et de l’innovation en statistique (DRIS), Statistique Canada.

Rubin-Bleuer, S., Yung, W. et Landry, S. (2010). Adjusted maximum likelihood method for a small area model accounting for time and area effects. SRID-2010-006E, Division de la recherche et de l’innovation en statistique (DRIS), Statistique Canada.

Rubin-Bleuer, S., Yung, W. et Landry, S. (2011). Adjusted maximum likelihood method for a small area model accounting for time and area effects. Long abstract, Small Area Estimation, (SAE 20122) à Trèves, Allemagne, International Statistical Institute Satellite Conference.

Rubin-Bleuer, S., Yung, W. et Landry, S. (2012). Variance Component Estimation through the Adjusted Maximum Likelihood Approach. Exposé donné à la conférence en l’honneur du 75^e anniversaire de J.N.K. Rao, Université Carleton, mai 2012, Ottawa.

Yoshimori, M., et Lahiri, P. (2014). A new adjusted maximum likelihood method for the Fay-Herriot small area model. Journal of Multivariate Analysis, 124, 281-294.

Yuan, P. (2009). Comparison of SAE methods of variance estimation. Document interne, Division de la recherche et de l’innovation en statistique (DRIS), Statistique Canada.

Yuan, K.H., et Jennrich, R. (1998). Asymptotics of estimating equations under natural conditions. Journal of Multivariate Analysis, 65, 2, 245-260.

• Précédent
• Suivant