Inférence bayésienne pour un modèle des composantes de la variance fondée sur la vraisemblance composite par paire à partir des données d’enquête
Section 3. Études par simulations

• Précédent
• Suivant

3.1 Plan de simulation

Au moyen d’études par simulations, nous avons évalué le rendement de la méthode proposée, c’est-à-dire la vraisemblance composite par paire avec un ajustement de la courbure, et nous l’avons comparé à celui de la vraisemblance complète et de la vraisemblance composite par paire. Nous avons utilisé le modèle de l’équation (1.1) pour générer nos données, c’est-à-dire que pour $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}n$ et $j\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,m,$ nous avons simulé les valeurs de $Y_{ij}$ à partir de

$$Y_{ij}\mathrm{<mo>=</mo>}\theta +u_i+e_{ij}\mathrm{,}(3.1)$$

où $\theta \mathrm{<mo>=</mo>}1,$ $u_i\stackrel{\text{iid}}{{\displaystyle ~}}N\left(\mathrm{0,}{\sigma }_u^2\right)$ et $e_{ij}\stackrel{\text{iid}}{{\displaystyle ~}}N\left(\mathrm{0,}{\sigma }_e^2\right)\mathrm{.}$ Cela équivaut à avoir appliqué la génération de la superpopulation et l’échantillonnage décrits au paragraphe se rapportant à l’équation (1.1).

Notre première étude, qui ne fait pas l’objet du présent article, s’est penchée sur l’inférence au sujet de $\theta$ avec ${\sigma }_u$ connue et ${\sigma }_e.$ Elle a démontré que l’utilisation de la vraisemblance composite par paire pour l’inférence au sujet de $\theta$ surestimait grandement la précision et que l’ajustement de la courbure était efficace. Nous avons donc procédé à une étude plus approfondie visant à examiner l’inférence à la fois pour $\theta$ et ${\sigma }_u.$ Par souci de simplicité, nous avons pris ${\sigma }_e\mathrm{<mo>=</mo>}\text{0,5}$ et avons considéré $n\in \left\{20;40\right\}$ et $m\in \left\{5;10\right\}$. Pour la loi a priori demi-Cauchy définie dans l’équation (2.4), nous avons pris $A\in \left\{5;10;15\right\}.$ Il y avait 500 ensembles de données qui se répliquaient pour chaque paramètre.

Nous avons examiné trois scénarios : 1) ${\sigma }_u\in \left\{\text{0,1;}\text{0,5}\right\}$ et la loi a priori demi-Cauchy sur ${\sigma }_u;$ 2) rapport signal-bruit, $\text{RSB}\in \left\{\text{0,25;}\text{0,75}\right\}$ et la loi a priori demi-Cauchy sur ${\sigma }_u,$ où le $\text{RSB}\mathrm{<mo>=</mo>}{\sigma }_u^2/\left({\sigma }_u^2+{\sigma }_e^2\right);$ et 3) ${\sigma }_u\in \left\{\text{0,1;}\text{0,5}\right\}$ et une loi a priori uniforme ${\sigma }_u.$ Pour l’ensemble des scénarios, nous avons utilisé une loi a priori uniforme pour $\theta .$

À la section 3.2, nous décrivons les algorithmes utilisés dans les études par simulations.

3.2  Algorithmes

Comme nous l’avons fait aux sections 2.1 et 2.2, définissons $y\left(n\right)\mathrm{<mo>=</mo>}\left\{y_1\mathrm{,}\dots \mathrm{,}{y}_n\right\}$ avec $y_i\mathrm{<mo>=</mo>}{\left(y_{i1}\mathrm{,}\dots \mathrm{,}{y}_{im}\right)}^{\text{T}}$ et $\overline{y}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^m{y}_{ij}/\left(mn\right)}}.$ De plus, ${\eta }^{\left(t\right)}$ désigne la valeur de $\eta$ à la $t^{\text{e}}$ itération où $\eta \mathrm{<mo>=</mo>}{\left(\theta \mathrm{,}{\sigma }_u\right)}^{\text{T}}.$ La vraisemblance complète est

$$L_{\text{VC}}\left(\theta \mathrm{,}{\sigma }_u\mathrm{<mo>|</mo>}y\left(n\right)\right)\propto {\left|{\Sigma }_m\right|}^{-n/2}\text{exp}\left[-\frac{1}{2}\text{tr}\left({\Sigma }_m^{-1}{S}_0\right)\right]\mathrm{,}(3.2)$$

comme dans l’équation (2.7).

L’utilisation de l’équation (3.2) de même que de la loi a priori, $\pi \left(\eta \right),$ permet d’obtenir la densité a posteriori,

$$p_{\text{VC}}\left(\eta |y\left(n\right)\right)\propto L_{\text{VC}}\left(\eta |y\left(n\right)\right)\pi \left(\eta \right)\mathrm{.}$$

L’échantillonnage de $\theta$ et de ${\sigma }_u$ est réalisé en deux étapes :

1. Étape 1. Échantillonner ${\theta }^{\left(t\right)}$ à partir de $p_{\text{VC}}\left(\theta |y\left(n\right)\mathrm{,}{\sigma }_u^{(t-1)}\right)$ où $$\theta |\left(y\left(n\right)\mathrm{,}{\sigma }_u\right)~N\left(\overline{y},\frac{{\sigma }_e^2+m{\sigma }_u^2}{mn}\right)\mathrm{.}$$ Nous établissons la valeur de départ, ${\sigma }_u^{\left(0\right)},$ comme l’estimation du maximum de vraisemblance de ${\sigma }_u.$
2. Étape 2. Utiliser l’algorithme de Metropolis-Hastings (MH) pour tirer ${\sigma }_u^{\left(t\right)}$ à partir de $p_{\text{VC}}\left({\sigma }_u|y\left(n\right)\mathrm{,}{\theta }^{(t)}\right).$ Cette dernière valeur est facilement obtenue à partir de $p_{\text{VC}}\left(\eta |y\left(n\right)\right).$ Soit $s\mathrm{<mo>></mo>}0,$ le candidat ${\sigma }_u,$ désigné ${\sigma }_u^{*},$ est tiré de la distribution des sauts, $N\left({\sigma }_u^{\left(t-1\right)}\mathrm{,}{s}^2\right).$ Si ${\sigma }_u^{*}<\mathrm{0,}{\sigma }_u^{\left(t\right)}\mathrm{<mo>=</mo>}{\sigma }_u^{\left(t-1\right)}.$ Sinon, la procédure est normale et comporte un rapport d’acceptation-rejet $p_{\text{VC}}\left({\eta }_{\text{VC}}^{*}|y\left(n\right)\right)/p_{\text{VC}}\left({\eta }_{\text{VC}}^{(t-1)}|y\left(n\right)\right)$ où $\eta {}_{\text{VC}}^{*}={\left({\theta }^{\left(t\right)}\mathrm{,}{\sigma }_u^{*}\right)}^{\text{T}}$ et $\eta {}_{\text{VC}}^{(t-1)}={\left({\theta }^{(t)}\mathrm{,}{\sigma }_u^{(t-1)}\right)}^{\text{T}}.$
3. Étape 3. Répéter les étapes 1 et 2 pour $K\mathrm{<mo>=</mo>}1000$ fois avec les 200 premiers échantillons utilisés comme rodage.

La vraisemblance composite par paire (VP) est

$$L_{\text{VP}}\left(\theta \mathrm{,}{\sigma }_u|y\left(n\right)\right)\propto {\left|\Sigma {}_2\right|}^{nm\left(m-1\right)/4}\text{exp}\left[-\frac{1}{2}\text{tr}\left({\Sigma }_2^{-1}{S}_{0\text{VP}}\right)\right]\mathrm{,}(3.3)$$

comme dans l’équation (2.8).

L’utilisation de l’équation (3.3) de même que de la loi a priori choisie, $\pi \left(\eta \right),$ permet d’obtenir la densité a posteriori $p_{\text{VP}}\left(\eta |y\left(n\right)\right).$

L’échantillonnage de $\theta$ et de ${\sigma }_u$ est réalisé en trois étapes :

1. Étape 1. Échantillonner ${\theta }^{\left(t\right)}$ à partir de $p_{\text{VP}}\left(\theta |y\left(n\right)\mathrm{,}{\sigma }_u^{(t-1)}\right)$ où $$\theta |\left(y\left(n\right)\mathrm{,}{\sigma }_u\right)~N\left(\overline{y},\frac{{\sigma }_e^2+2{\sigma }_u^2}{nm\left(m-1\right)}\right)\mathrm{.}$$
2. Étape 2. Utiliser l’algorithme de Metropolis-Hastings (MH) pour tirer ${\sigma }_u^{\left(t\right)}$ à partir de $p_{\text{VP}}\left({\sigma }_u|y\left(n\right)\mathrm{,}{\theta }^{(t)}\right),$ tel qu’il est décrit à l’étape 2 ci-dessus pour la VC (en remplaçant la VP par la VC dans toutes les formules).
3. Étape 3. Répéter les étapes 1 et 2 pour $K\mathrm{<mo>=</mo>}1000$ fois avec les 200 premiers échantillons utilisés comme rodage.

La dernière étape consiste à obtenir la vraisemblance composite par paire ajustée (VPA) (courbure), tel qu’elle est décrit à la section 2.3. Cette dérivation, fondée sur l’approche de Ribatet, Cooley et Davison, exploite ${\widehat{\eta }}_{\text{VPA}},$ les moyennes a posteriori estimées de $\theta$ et de ${\sigma }_u.$

1. Étape 1. Soit $\left(s_{\theta }\mathrm{,}{s}_{\sigma }\right),$ tirer le candidat ${\eta }^{*}\mathrm{<mo>=</mo>}{\left({\theta }^{*}\mathrm{,}{\sigma }_u^{*}\right)}^{\text{T}}$ à partir de la distribution des sauts normale à deux variables, $N_2\left(\eta {}^{\left(t-1\right)}\mathrm{,}\Sigma \right)$ où $\Sigma \mathrm{<mo>=</mo>}\text{diag}\left(s_{\theta }^2\mathrm{,}{s}_{\sigma }^2\right).$ Si ${\sigma }_u^{*}\mathrm{<mo><</mo>0,}{\eta }^{\left(t\right)}\mathrm{<mo>=</mo>}{\eta }^{\left(t-1\right)}.$ Sinon, passer à l’étape 2.
2. Étape 2. Définissons $l_{\text{VP}}\left(y\left(n\right)|\theta \mathrm{,}{\sigma }_u\right)$ comme le logarithme du rapport de vraisemblance composite par paire obtenu en prenant le logarithme de l’équation (3.3) et $l_{\text{VP}}\left(y_i|\theta \mathrm{,}{\sigma }_u\right)$ comme le logarithme du rapport de vraisemblance composite par paire correspondant aux données de la grappe $i,$ c’est-à-dire $y_i.$
3. Étape 3. Obtenir numériquement $\widehat{H}\mathrm{<mo>=</mo>}{\nabla }^2{l}_{\text{VP}}\left(y\left(n\right)|{\widehat{\theta }}_{\text{VP}},{\widehat{\sigma }}_{u\text{VP}}\right)$ et $$\widehat{J}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n\left[\nabla l_{\text{VP}}\left(y_i|{\widehat{\theta }}_{\text{VP}}\mathrm{,}{\widehat{\sigma }}_{u\text{VP}}\right){\left\{\nabla l\left(y_i|{\widehat{\theta }}_{\text{VP}}\mathrm{,}{\widehat{\sigma }}_{u\text{VP}}\right)\right\}}^{\text{T}}\right]}\mathrm{,}$$ où ${\widehat{\theta }}_{\text{VP}}$ et ${\widehat{\sigma }}_{u\text{VP}}$ sont les moyennes a posteriori estimées $\theta$ et ${\sigma }_u.$
4. Étape 4. D’après l’approche de Ribatet, Cooley et Davison, et au moyen de la décomposition en valeurs singulières, nous écrivons $\widehat{H}\mathrm{<mo>=</mo>}{M}^{\text{T}}M$ et $\widehat{H}{\widehat{J}}^{-1}\widehat{H}\mathrm{<mo>=</mo>}{M}_A^{\text{T}}{M}_A$ pour certaines matrices $M$ et $M_A.$ Définir ensuite $C\mathrm{<mo>=</mo>}{M}^{-1}{M}_A.$ Dans notre cas, $C$ est une matrice $2\times 2.$
5. Étape 5. Selon l’approche de Ribatet, Cooley et Davison, le logarithme du rapport de vraisemblance composite par paire ajusté, $l_{\text{VPA}},$ est $$l_{\text{VPA}}\left(y\left(n\right)|\eta \right)\mathrm{<mo>=</mo>}{l}_{\text{VP}}\left(y\left(n\right)|{\eta }^{*}\right)$$ où $${\eta }^{*}\mathrm{<mo>=</mo>}{\widehat{\eta }}_{\text{VP}}+C\left(\eta -{\widehat{\eta }}_{\text{VP}}\right).$$
6. Étape 6. Définir la densité a posteriori par paire ajustée comme $$p_{\text{VPA}}\left(\eta |y\left(n\right)\right)\propto L_{\text{VPA}}\left(y\left(n\right)|\eta \right)\pi \left(\theta \mathrm{,}{\sigma }_u\right)$$ où $L_{\text{VPA}}\left(y\left(n\right)|\eta \right)\mathrm{<mo>=</mo>}\text{exp}\left(l_{\text{VPA}}\left(y\left(n\right)|\eta \right)\right),$ ce dernier étant défini à l’étape 5.
   À l’aide de la valeur candidate, ${\eta }^{*},$ obtenue à l’étape 1, définir la valeur candidate ajustée ${\eta }_c^{*}\mathrm{<mo>=</mo>}{\widehat{\eta }}_{\text{VP}}+C\left({\eta }^{*}-{\widehat{\eta }}_{\text{VP}}\right).$ Le rapport d’acceptation-rejet est alors $$p_{\text{VPA}}\left({\text{η}}_c^{*}|y\left(n\right)\right)/p_{\text{VPA}}\left({\eta }^{\left(t\right)}|y\left(n\right)\right).$$ Les autres étapes sont les étapes normales pour l’algorithme de Metropolis-Hastings.

3.3 Résultats des simulations

Pour chaque méthode (VC, VP, VPA), chaque paramètre du plan $\left(m\mathrm{,}n\right)$ et chaque loi a priori, nous résumons les résultats des simulations en utilisant a) les taux de couverture des intervalles de crédibilité sous échantillonnage répété et b) les moyennes des points 0,025; 0,25; 0,50; 0,75 et 0,975 des lois a posteriori de $\theta$ et de ${\sigma }_u.$

Des résumés graphiques sont également présentés, c’est-à-dire des estimations de la densité a posteriori moyenne pour chacune des lois a posteriori, c’est-à-dire $p_{\text{VC}}\left(\eta |y\left(n\right)\right),p_{\text{VP}}\left(\eta |y\left(n\right)\right)$ et $p_{\text{VPA}}\left(\eta |y\left(n\right)\right).$ D’abord, considérons un intervalle, disons, $\left[a\mathrm{,}b\right],$ qui soutient la plus grande partie de la masse (par exemple 95 %) des densités a posteriori. Divisons-le ensuite en $M\mathrm{<mo>=</mo>}50$ sous-intervalles échelonnés également comportant les points de découpage $a\mathrm{<mo>=</mo>}{c}_0<c_1<\dots <c_{M-1}<c_M\mathrm{<mo>=</mo>}b.$ Pour $t\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,T,$ soit ${\widehat{f}}_P^{(t)}\left(\mathrm{.}\right)$ l’estimation de la densité a posteriori $f_P\left(\mathrm{.}\right),$ dérivée de la simulation $t^{\text{e}},$ où $P$ désigne la VC, la VP ou la VPA et $T$ correspond au nombre de simulations. Définissons ensuite, pour $r\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,M,$

$${\widehat{f}}_P\left(c_r\right)\mathrm{<mo>=</mo>}\frac{1}{T}{\displaystyle \sum _{t\mathrm{<mo>=</mo>}1}^T{\widehat{f}}_P^{\left(t\right)}\left(c_r\right)}\mathrm{.}$$

Prenons ensuite une courbe reliant les points $\left\{c_r\mathrm{,}{\widehat{f}}_P\left(c_r\right)\right\}$ pour $a\mathrm{<mo>=</mo>}{c}_0<c_1<\dots <c_M\mathrm{<mo>=</mo>}b$ comme estimation de la densité a posteriori moyenne pour $f_P\left(\mathrm{.}\right).$

Le tableau 3.1 présente les taux de couverture pour $\theta$ et ${\sigma }_u$ pour $A\mathrm{<mo>=</mo>}\mathrm{15,}$ $n\in \left\{20;40\right\},$ $m\in \left\{5;10\right\}$ et ${\sigma }_u\in \left\{\text{0,1;}\text{0,289;}\text{0,5;}\text{0,866}\right\}\mathrm{.}$ La figure 3.1 présente les estimations moyennes de la densité a posteriori pour $\theta$ et ${\sigma }_u$ pour $A\mathrm{<mo>=</mo>}\mathrm{15,}{\sigma }_u\in \left\{\text{0,1;}\text{0,5}\right\}\mathrm{,}n\mathrm{<mo>=</mo>}40$ et $m\mathrm{<mo>=</mo>}10.$ Au tableau 3.1 et à la figure 3.1, les résumés sont présentés pour la vraisemblance complète (VC), la vraisemblance composite par paire (VP) et la vraisemblance composite par paire ajustée (VPA).

Le résumé suivant présente les résultats pour la loi a priori demi-Cauchy seulement avec $A\in \left\{5;10;15\right\},$ $m\in \left\{5;10\right\},$ $n\in \left\{20;40\right\}$ et ${\sigma }_u\in \left\{\text{0,1;}\text{0,289;}\text{0,5;}\text{0,866}\right\},$ les deuxième et quatrième valeurs correspondant au RSB = 0,25 et au RSB = 0,75, respectivement. Les résultats sont semblables pour les trois choix de $A$ et pour la loi a priori uniforme.

Sans ajustement, les couvertures de la VP sont considérablement différentes du niveau nominal de 0,95. Par exemple (tableau 3.1), pour $A\mathrm{<mo>=</mo>}15,$ $n\mathrm{<mo>=</mo>}40,$ $m\mathrm{<mo>=</mo>}10$ et ${\sigma }_u\mathrm{<mo>=</mo>}\text{0,5},$ la couverture pour $\theta$ est inférieure à 0,30. Si nous considérons toutes les valeurs des paramètres du plan, la plus grande couverture est de 0,70. Dans la plupart des cas, la couverture pour $\theta$ est bien inférieure à 0,70.

Lorsque l’ajustement de la courbure est pris en compte, la couverture pour $\theta$ est excellente. Des 48 cas (trois choix de $A,$ deux choix de $m,$ deux choix de $n,$ quatre choix de ${\sigma }_u),$ treize avaient une couverture entre 0,93 et 0,95, vingt-deux entre 0,92 et 0,93, onze entre 0,91 et 0,92, et deux de moins de 0,91, avec le dernier cas de figure pour ${\sigma }_u\mathrm{<mo>=</mo>}\text{0,1},$ $n\mathrm{<mo>=</mo>}20,$ $m\mathrm{<mo>=</mo>}10$ et $A=5$ et 15.

Lorsque l’ajustement de la courbure est pris en compte, la couverture pour ${\sigma }_u$ varie considérablement, mais il y a, dans la plupart des cas, une très grande amélioration de la couverture comparativement à l’utilisation de la vraisemblance composite par paire non ajustée.

Les tracés (figure 3.1) montrent que, pour $\theta ,$ la loi a posteriori correspondant à la vraisemblance ajustée se rapproche beaucoup de la loi a posteriori fondée sur la vraisemblance complète. Pour ${\sigma }_u,$ il existe des différences entre les lois a posteriori correspondant à la vraisemblance complète et à la vraisemblance ajustée, en particulier une transition vers des valeurs plus faibles dans le cas de cette dernière.

Pour étudier les effets de l’augmentation de $m$ et de $n,$ considérons la différence $\delta \mathrm{<mo>=</mo>}{C}_{\text{VC}}-C_{\text{VPA}},$ où $C$ désigne la couverture et VC et VPA renvoient aux lois a posteriori correspondantes.

Dans l’ensemble, lorsque nous considérons tous les $m\mathrm{,}n\mathrm{,}A$ et ${\sigma }_u,$ pour $\theta ,$ $\delta$ diminue à mesure que $n$ augmente. Pour les valeurs plus grandes de ${\sigma }_u,$ $\delta$ diminue à mesure que $m$ augmente, tandis que pour les valeurs plus faibles de ${\sigma }_u,$ $\delta$ a tendance à augmenter à mesure que $m$ augmente. Dans l’ensemble, pour ${\sigma }_u,$ $\delta$ diminue à mesure que $n$ augmente, sauf dans le cas où ${\sigma }_u\mathrm{<mo>=</mo>}\text{0,1},$ tandis que $\delta$ augmente à mesure que $m$ augmente.

La raison de la détérioration de l’ajustement à mesure que $m$ augmente pourrait être que le nombre de paires par grappe est $m\left(m-1\right)/2$ et augmente plus rapidement, de sorte que la vraisemblance par paire se concentre plus rapidement autour de son mode; l’ajustement de la courbure peut ne pas être suffisant pour compenser un changement de forme du logarithme du rapport de vraisemblance composite par paire, par exemple une augmentation de l’aplatissement.

Le tableau 3.2 présente les taux de non-couverture unilatéraux pour les 95 % des intervalles de crédibilité pour $\theta$ et ${\sigma }_u$ au moyen de $A\mathrm{<mo>=</mo>}15,$ $n\in \left\{20;40\right\},$ $m\in \left\{5;10\right\}$ et ${\sigma }_u\in \left\{\text{0,1;}\text{0,289;}\text{0,5;}\text{0,866}\right\}.$ Nous observons ce qui suit :

1. Pour $\theta ,$ la non-couverture pour les intervalles de vraisemblance complète semble symétrique. La vraisemblance par paire ajustée présente un sous-dénombrement pour $\theta$ et, sauf lorsque ${\sigma }_u$ est égale à 0,1, le sous-dénombrement est symétrique. Nous observons seulement une dépendance de la couverture sur $m$ dans le cas où ${\sigma }_u\mathrm{<mo>=</mo>}\text{0,1}.$
2. Pour ${\sigma }_u,$ l’intervalle de vraisemblance complète présente un sous-dénombrement qui se rapproche du niveau nominal et n’est pas très asymétrique, sauf dans le cas où ${\sigma }_u\mathrm{<mo>=</mo>}\text{0,1},$ pour lequel nous observons un sous-dénombrement marqué. Pour ${\sigma }_u>\text{0,1}$ et $m\mathrm{<mo>=</mo>}5,$ la couverture s’améliore à mesure que $n$ passe de 20 à 40, mais pour ${\sigma }_u>\text{0,1}$ et $m\mathrm{<mo>=</mo>}10,$ la couverture présente peu de différence pour les deux valeurs de $n.$
3. Pour ${\sigma }_u,$ la vraisemblance par paire ajustée présente une couverture asymétrique. Sauf dans le cas où ${\sigma }_u\mathrm{<mo>=</mo>}\text{0,1},$ l’ampleur de la non-couverture à gauche a tendance à être semblable à celle de la vraisemblance complète, mais beaucoup plus grande à droite, et la couverture s’améliore à mesure que $n$ passe de 20 à 40.

Sachant que le logarithme du rapport de vraisemblance par paire ajusté n’est pas explicitement construit pour estimer le logarithme du rapport de vraisemblance complet, nous constatons que, dans la figure 3.1, le logarithme du rapport de vraisemblance par paire ajusté diminue plus rapidement dans les queues.

Nous avons également tenté de centrer l’ajustement de la courbure sur le mode a posteriori du logarithme du rapport de vraisemblance par paire plutôt que sur la moyenne a posteriori du logarithme du rapport de vraisemblance par paire, et nous avons constaté que le sous-dénombrement augmentait, même si l’asymétrie de la couverture était moins marquée, pour les intervalles de crédibilité qui en découlent.

• Précédent
• Page principale de l'article
• Suivant