La vraisemblance pénalisée de Firth pour les régressions à risques proportionnels en cas d’enquêtes complexes
Section 1. Introduction

• Précédent
• Suivant

Le modèle de régression à risques proportionnels de Cox (Cox, 1972) est couramment utilisé dans l’analyse des données de survie. Il s’agit d’un modèle semi-paramétrique qui explique l’effet des variables explicatives sur les taux de risque. Le modèle suppose que l’effet des variables explicatives a une forme linéaire, mais il permet que la fonction de survie sous-jacente ait une forme non spécifiée. On estime les paramètres du modèle en maximisant une vraisemblance partielle (Cox, 1972, 1975).

Pour estimer les paramètres canoniques dans les distributions de la famille exponentielle, Firth (1993) a proposé de multiplier la vraisemblance par la loi a priori de Jeffreys afin d’obtenir une estimation par le maximum de vraisemblance qui soit du premier ordre sans biais. La vraisemblance pénalisée prend la forme

$$L_p\left(\beta \right)\mathrm{<mo>=</mo>}L\left(\beta \right){\left|I\left(\beta \right)\right|}^{\text{0,5}}$$

où $L\left(\beta \right)$ est la vraisemblance non pénalisée, $I$ est la matrice d’information et $\beta$ est un vecteur des paramètres de régression. La vraisemblance pénalisée de Firth est une technique très utile en pratique, non seulement aux fins de réduction du biais, mais aussi de correction des vraisemblances monotones.

Souvent, les modèles de régression à risques proportionnels souffrent de vraisemblances monotones, dans lesquelles la vraisemblance converge vers une valeur finie, mais un paramètre au moins diverge (Heinze, 1999). La vraisemblance pénalisée de Firth sert également à corriger les vraisemblances monotones et à obtenir des estimations de paramètres qui convergent (Heinze, 1999; Heinze et Schemper, 2001; Heinzel, Rüdiger et Schilling, 2002).

Bien que la vraisemblance pénalisée de Firth soit utile aux fins de réduction des biais et d’obtention d’estimations à partir de vraisemblances monotones, elle n’a pas été étudiée sur des enquêtes complexes comportant des poids inégaux. Il est raisonnable d’utiliser une vraisemblance pondérée pour les enquêtes complexes afin de compenser la pondération inégale (Fuller, 1975; Binder et Patak, 1994). Les ensembles de données d’enquête comprennent couramment des poids de sondage ou des poids d’analyse pour lesquels la somme des poids est un estimateur de la taille de la population. Toutefois, ces pondérations non mises à l’échelle ne mettent pas adéquatement à l’échelle la matrice d’information utilisée dans le terme de pénalité. Il est souhaitable que les paramètres d’une régression à risques proportionnels pour données d’enquête aient les deux propriétés suivantes.

• Invariance : les estimations ponctuelles et les erreurs-types pour les paramètres de régression doivent être invariantes par rapport à l’échelle des poids.
• Précision : la variance linéarisée de Taylor pour les paramètres de régression estimés doit être proche de la variance du jackknife avec suppression d’une UPE.

Dans l’article, nous montrons d’abord que si l’on n’utilise pas la correction de Firth, alors l’invariance et la précision sont satisfaites, mais que si la correction de Firth est utilisée avec les poids non mis à l’échelle, alors les estimations ponctuelles et les erreurs-types ne sont pas invariantes par rapport à l’échelle des poids. Autrement dit, si les poids sont multipliés par une constante et que la correction de Firth est utilisée, les estimations ponctuelles et les erreurs-types seront différentes. Nous proposons ensuite une méthode de mise à l’échelle des poids de sens commun pour démontrer que la correction de Firth utilisant des poids mis à l’échelle possède les deux propriétés souhaitables. La seule différence entre les poids mis à l’échelle et ceux non mis à l’échelle est que la somme des poids mis à l’échelle est égale à la taille de l’échantillon, alors que la somme des poids non mis à l’échelle est un estimateur de la taille de la population.

1.1  Exemple d’utilisation de poids non mis à l’échelle

Nous avons utilisé un ensemble de données tiré d’une étude portant sur 65 patients atteints de myélome qui ont été traités au moyen d’agents alkylants (Lee, Wei et Amato, 1992) pour démontrer les propriétés de la vraisemblance pénalisée de Firth utilisant des poids non mis à l’échelle. Les durées de survie en mois ont été enregistrées pour chaque patient. Les patients qui étaient vivants après la période de l’étude étaient considérés comme des données censurées. Nous disposions des variables suivantes pour chaque patient :

• Durée [Time] : durée de la survie en mois,
• Statut vital [Vstatus] : état du patient, zéro ou un, indiquant respectivement si le patient était vivant ou mort,
• LogBUN : log du niveau d’azote uréique sanguin,
• HGB : taux d’hémoglobine dans le sang.

Pour créer une vraisemblance monotone, nous avons ajouté une nouvelle variable explicative, Contrived [Artificielle], telle que sa valeur à tout moment de l’événement est la plus grande de toutes les valeurs de l’ensemble de risques (voir l’exemple « Correction de Firth pour vraisemblance monotone » dans « The PHREG Procedure » [La procédure PHREG] dans SAS Institute Inc. (2018)). La variable Contrived [Artificielle] a la valeur 1 si la durée de survie observée est inférieure ou égale à 65; sinon, elle a la valeur 0.

Pour démontrer l’effet des poids dans la vraisemblance pénalisée de Firth, nous avons créé trois variables de poids, $w1,$ $w3$ et $w5,$ avec respectivement des valeurs de 1, 1 000 et 100 000 pour chaque observation. On estime les paramètres de régression à risques proportionnels en maximisant une vraisemblance pondérée comme cela est décrit dans la section 1.2. Parce que $w1$ a la valeur 1 pour toutes les observations, l’utilisation de $w1$ dans l’analyse équivaut à l’exécution de l’analyse non pondérée.

Nous avons ajusté les deux modèles à risques proportionnels suivants à l’aide de la procédure PHREG dans SAS/STAT^MD (voir « The PHREG Procedure » dans SAS Institute Inc. (2018) :

$$\lambda \left(t\mathrm{,}Z\right)\mathrm{<mo>=</mo>}{\lambda }_0\left(t\right)\exp \left({\beta }_1\text{LogBUN}+{\beta }_2\text{HGB}\right)$$

$$\lambda \left(t\mathrm{,}Z\right)\mathrm{<mo>=</mo>}{\lambda }_0\left(t\right)\exp \left({\beta }_1\text{LogBUN}+{\beta }_2\text{HGB}+{\beta }_3\text{Contrived}\right)$$

où $\lambda \left(t\right)$ et ${\lambda }_0\left(t\right)$ sont respectivement la fonction de risque et la fonction de risque de référence. La vraisemblance pénalisée de Firth n’est pas requise pour l’ajustement du premier modèle sans la variable Contrived (la probabilité converge en trois étapes d’itération), mais le deuxième modèle contenant la variable Contrived ne converge pas sans la pénalité de Firth dans la vraisemblance. Le tableau 1.1 présente la valeur de la vraisemblance et les trois coefficients de régression pour 14 itérations. Bien que la fonction objective et les coefficients de LogBUN et HGB convergent vers une valeur finie après la quatrième itération, ceux de Contrived divergent. Il s’agit d’un exemple de vraisemblance monotone pour la variable Contrived. En raison de cette monotonie, il faut utiliser la vraisemblance pénalisée de Firth pour ajuster le deuxième modèle contenant Contrived.

Si Contrived n’est pas utilisée comme variable explicative, les trois ensembles de poids produisent des estimations ponctuelles et des estimations de variance linéarisées de Taylor identiques (tableau 1.2). Les estimations de la variance du jackknife avec suppression d’une UPE sont également identiques pour les trois ensembles de poids. Ainsi, les estimations ponctuelles et les erreurs-types sont invariantes par rapport à l’échelle des poids quand la correction de Firth n’est pas utilisée.

En revanche, si l’on utilise des poids non mis à l’échelle, les estimations ponctuelles pour Contrived ne sont pas invariantes par rapport à l’échelle des poids. Le tableau 1.3 présente les estimations des paramètres pour trois ensembles de poids quand Contrived est utilisée comme variable explicative (et que la vraisemblance pénalisée de Firth est appliquée). Parce que la vraisemblance n’est pas monotone (tableau 1.1) pour LogBUN et HGB, les estimations ponctuelles pour ces deux coefficients ne sont pas affectées par l’échelle des poids.

Si Contrived n’est pas utilisée comme variable explicative, le rapport entre les erreurs-types du jackknife et les erreurs-types de la linéarisation de Taylor est de 1,13 et 1,10 pour les trois ensembles de poids des variables LogBUN et HGB, respectivement. Ainsi, le rapport entre la variance par la méthode du jackknife et la variance linéarisée de Taylor pour la vraisemblance non pénalisée est invariant par rapport à l’échelle des poids, et il est raisonnable de penser que le rapport est invariant quand on utilise la vraisemblance pénalisée.

1.2  Bref examen des estimations ponctuelles et des estimations de variance pour les paramètres de régression de populations finies

Avant de discuter de la méthode de mise à l’échelle des poids, nous examinerons brièvement les estimations ponctuelles et les estimations de variance pour les paramètres de régression dans une régression à risques proportionnels sur des enquêtes complexes comportant des poids inégaux. Lin et Wei (1989), Binder (1990, 1992), Lin (2000) et Boudreau et Lawless (2006) ont traité de l’estimation par la méthode du pseudo-maximum de vraisemblance des paramètres de régression à risques proportionnels pour les données d’enquête. On trouve une description plus générale de l’estimation des paramètres de régression pour les enquêtes complexes dans Kish et Frankel (1974); Godambe et Thompson (1986); Pfeffermann (1993), Korn et Graubard (1999, chapitre 3), Chambers et Skinner (2003, chapitre 2) et Fuller (2009, section 6.5). Wolter (2007) a décrit plusieurs techniques d’estimation de la variance pour les données d’enquête.

Soit $U_N\mathrm{<mo>=</mo>}\left\{\mathrm{1,}\mathrm{2,}\dots \mathrm{,}N\right\}$ l’ensemble des indices et $F_N$ l’ensemble des valeurs pour une population finie de taille $N.$ On suppose que la durée de survie de chaque membre de la population finie suit sa propre fonction de risque, ${\lambda }_i\left(t\right),$ exprimée comme suit :

$${\lambda }_i\left(t\right)\mathrm{<mo>=</mo>}\lambda \left(t\mathrm{<mo>;</mo>}{Z}_i\left(t\right)\right)\mathrm{<mo>=</mo>}{\lambda }_0\left(t\right)\exp \left(Z_i^{\prime }\left(t\right)\beta \right)$$

où ${\lambda }_0\left(t\right)$ est une fonction de risque de référence arbitraire et non spécifiée, $Z_i\left(t\right)$ est un vecteur de taille $P$ des variables explicatives de l’unité $i^{\text{e}}$ au temps $t,$ et $\beta$ est un vecteur de paramètres de régression inconnus.

La fonction de vraisemblance partielle introduite par Cox (1972, 1975) élimine le risque de référence inconnu ${\lambda }_0\left(t\right)$ et tient compte des durées de survie censurées. Si toute la population est observée, cette fonction de vraisemblance partielle peut servir à estimer $\beta .$ Soit ${\beta }_N$ l’estimateur souhaité.

En supposant un modèle de travail avec des réponses non corrélées, on obtient ${\beta }_N$ en maximisant la log-vraisemblance partielle,

$$l_N\left(\beta \right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\in U_N}}\log \left\{\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)\right\}$$

par rapport à $\beta ,$ où $\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)$ est la fonction de vraisemblance partielle de Cox.

Supposons qu’un échantillon probabiliste $A_N$ est sélectionné dans la population finie $U_N.$ Soit ${\pi }_i$ la probabilité de sélection et $w_i\left(={\pi }_i^{-1}\right)$ le poids d’échantillonnage pour l’unité $i.$ Supposons ensuite que les variables explicatives $Z_i\left(t\right)$ et la durée de survie $t_i$ sont disponibles pour chaque unité de l’échantillon $A_N.$ Un estimateur sans biais par rapport au plan pour la log-vraisemblance de la population finie est

$$l\left(\beta \right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\in A_N}}{\pi }_i^{-1}\log \left\{\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)\right\}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\in A_N}}{w}_i\log \left\{\text{L}\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)\right\}.$$

On peut obtenir un estimateur basé sur un échantillon ${\widehat{\beta }}_N$ pour la quantité finie de population ${\beta }_N$ en maximisant la pseudo-log-vraisemblance partielle $l\left(\beta \mathrm{<mo>;</mo>}{Z}_i\left(t\right)\mathrm{,}{t}_i\right)$ par rapport à $\beta .$ On obtient la variance fondée sur le plan pour ${\widehat{\beta }}_N$ en supposant que l’ensemble de valeurs de la population finie $F_N$ est fixe.

La vraisemblance de Breslow pondérée peut être exprimée comme suit :

$$L\left(\beta \right)\mathrm{<mo>=</mo>}{\displaystyle \prod _{k\mathrm{<mo>=</mo>\; 1}}^K}\frac{\exp \left({\beta }^{\prime }{\displaystyle {\sum }_{D_k}{w}_i}{Z}_i\left(t\right)\right)}{{\left\{{\displaystyle {\sum }_{R_k}{w}_i}\exp \left({\beta }^{\prime }{Z}_i\left(t\right)\right)\right\}}^{{\displaystyle {\sum }_{D_k}{w}_i}}}$$

où $R_k$ est l’ensemble de risques juste avant le $k^{\text{e}}$ temps d’événement ordonné $t_{\left(k\right)},$ $D_k$ est l’ensemble des personnes qui fait défaillance au temps $t_{\left(k\right)},$ et $K$ est le nombre de temps d’événements distincts.

On obtient les estimations ponctuelles pour $\beta$ en maximisant $l\left(\beta \right)\mathrm{<mo>=</mo>}\log \left[L\left(\beta \right)\right].$

Bien que les poids suffisent aux fins d’estimation des coefficients de régression pour la population finie, on doit aussi utiliser les données de stratification et de corrélation intra-grappe pour estimer la variabilité d’échantillonnage. Afin d’estimer la variabilité d’échantillonnage, on peut utiliser la méthode de linéarisation en séries de Taylor ou une méthode de rééchantillonnage.

1.2.1  Estimateur de la variance analytique par la méthode de linéarisation en séries de Taylor

La méthode de linéarisation en séries de Taylor utilise la somme des carrés des scores résiduels pondérés pour estimer la variabilité d’échantillonnage.

Définissons $\overline{Z}\left(\beta \mathrm{,}t\right)\mathrm{<mo>=</mo>}\frac{S^{\left(1\right)}\left(\beta \mathrm{,}t\right)}{S^{\left(0\right)}\left(\beta \mathrm{,}t\right)},$ où

$$S^{\left(0\right)}\left(\beta ,t\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{A_N}}{w}_{i}I\left(t_i\ge t\right)\exp \left({\beta }^{\prime }{Z}_i\left(t\right)\right)$$

et

$$S^{\left(1\right)}\left(\beta ,t\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{A_N}}{w}_{i}I\left(t_i\ge t\right)\exp \left({\beta }^{\prime }{Z}_i\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}\right)Z_i\left(t\right).$$

Le score résiduel pour le $i^{\text{e}}$ sujet est

$$\begin{array}{ll}{u}_i\left(\beta \right)=\hfill & {\Delta }_i\left\{Z_i\left(t_i\right)-\overline{Z}\left(\beta ,t_i\right)\right\}\hfill \\ \hfill & -{\displaystyle \sum _{j\in A_N}}\left[{\Delta }_j\frac{w_{j}I\left(t_i\ge t_j\right)\exp \left({\beta }^{\prime }{Z}_i\left(t_j\right)\right)}{S^{\left(0\right)}\left(\beta ,t_j\right)}\left\{Z_i\left(t_j\right)-\overline{Z}\left(\beta \mathrm{,}{t}_j\right)\right\}\right]\hfill \end{array}$$

où ${\Delta }_i$ est l’indicateur d’événement.

Alors, l’estimateur de variance linéarisé de Taylor est

$$\widehat{V}\left(\widehat{\beta }\right)\mathrm{<mo>=</mo>}{I}^{-1}\left(\widehat{\beta }\right)\text{G}{I}^{-1}\left(\widehat{\beta }\right)$$

où $I\left(\widehat{\beta }\right)$ est la matrice d’information observée et la $p\times p$ matrice $G$ est définie comme étant

$$G\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{,}j\in A_N\mathrm{<mo>:</mo>}i\mathrm{<mo><</mo>}j}}\frac{{\pi }_i{\pi }_j-{\pi }_{ij}}{{\pi }_{ij}}{\left(\frac{{\widehat{u}}_i}{{\pi }_i}-\frac{{\widehat{u}}_j}{{\pi }_j}\right)}^{\prime }\left(\frac{{\widehat{u}}_i}{{\pi }_i}-\frac{{\widehat{u}}_j}{{\pi }_j}\right)$$

où ${\pi }_{ij}$ sont les probabilités d’inclusion conjointe pour les unités $i$ et $j.$

En particulier, pour les plans d’échantillonnage en grappes stratifiés dans lesquels les UPE sont sélectionnées au moyen d’un échantillon aléatoire simple sans remise, la $p\times p$ matrice $G$ se réduit à

$$G\mathrm{<mo>=</mo>}{\displaystyle \sum _{h\mathrm{<mo>=</mo>\; 1}}^H}\frac{n_h\left(1-f_h\right)}{n_h-1}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^{n_h}}{\left(e_{hi+}-{\overline{e}}_{h\cdot \cdot }\right)}^{\prime }\left(e_{hi+}-{\overline{e}}_{h\cdot \cdot }\right)$$

où $e_{hi+}$ est la somme pondérée des scores résiduels, ${\widehat{u}}_{hij},$ dans la strate $h$ et l’UPE $i;$ ${\overline{e}}_{h\mathrm{..}}$ est la moyenne de $e_{hi+};$ $n_h$ est le nombre d’UPE; et $f_h$ est la fraction d’échantillonnage dans la strate $h.$

Ces estimateurs sont largement étudiés par la littérature sur les enquêtes-échantillons. Par exemple, Binder (1992) et Lin (2000) fournissent des conditions en vertu desquelles $\widehat{\beta }$ et $\widehat{V}\left(\widehat{\beta }\right)$ sont convergents. Chambless et Boyle (1985) ont calculé la variance fondée sur le plan et la normalité asymptotique pour des modèles à risques proportionnels discrets.

1.2.2  Estimateur de la variance par répliques au moyen de la méthode du jackknife avec suppression d’une UPE

La méthode du jackknife est une méthode d’estimation de la variance par répliques couramment utilisée en cas d’enquêtes complexes. Pour créer des rééchantillonnages, elle supprime (en lui attribuant un poids nul) une UPE à la fois de l’échantillon complet. Dans chaque rééchantillonnage, les poids d’échantillonnage des UPE restantes sont modifiés par le coefficient du jackknife ${\alpha }_r.$ Les poids modifiés sont appelés poids de rééchantillonnage.

Supposons que l’UPE $i_r$ dans la strate $h_r$ soit omise du $r^{\text{e}}$ rééchantillonnage, alors les poids de rééchantillonnage et les coefficients du jackknife sont donnés par

$$w_{hij}^{\left(r\right)}\mathrm{<mo>=</mo>}\{\begin{array}{ll}0\hfill & i\mathrm{<mo>=</mo>}{i}_r\text{et}h\mathrm{<mo>=</mo>}{h}_r\hfill \\ w_{hij}/{\alpha }_r\hfill & i\ne i_r\text{et}h\mathrm{<mo>=</mo>}{h}_r\hfill \\ w_{hij}\hfill & h\ne h_r\hfill \end{array}$$

et ${\alpha }_r\mathrm{<mo>=</mo>}\frac{n_{h_r}-1}{n_{h_r}},$ respectivement, pour toutes les unités d’observation $j$ dans la strate $h$ et l’UPE $i.$ Le nombre d’UPE dans la strate $h_r$ est de $n_{h_r}.$

On peut appliquer la méthode du jackknife pour estimer les variances des paramètres de régression estimés pour le modèle de Cox parce que les paramètres du modèle sont les solutions d’un ensemble d’équations d’estimation qui sont des fonctions lisses de totaux (les fonctions de score correspondantes sont données dans la section 2). Les propriétés des estimateurs de la variance par la méthode du jackknife pour les modèles de régression à risques proportionnels sont étudiées dans Shao et Tu (1995, section 8.3).

Pour appliquer la méthode du jackknife, on estime les paramètres du modèle au moyen de l’échantillon complet et en utilisant chaque échantillon répété. Soit $\widehat{\beta }$ les coefficients de régression à risques proportionnels estimés à partir de l’échantillon complet et soit ${\widehat{\beta }}_r$ les coefficients de régression estimés à partir du $r^{\text{e}}$ rééchantillonnage. Alors, la matrice de covariance de $\widehat{\beta }$ est estimée par

$$\widehat{V}\left(\widehat{\beta }\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R}{\alpha }_r\left({\widehat{\beta }}_r-\widehat{\beta }\right){\left({\widehat{\beta }}_r-\widehat{\beta }\right)}^{\prime }.$$

Si les fractions de sondage ne sont pas ignorables, la matrice de covariance de $\widehat{\beta }$ est estimée par

$$\widehat{V}\left(\widehat{\beta }\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>\; 1}}^R}{\alpha }_r\left(1-f_r\right)\left({\widehat{\beta }}_r-\widehat{\beta }\right){\left({\widehat{\beta }}_r-\widehat{\beta }\right)}^{\prime }$$

où $f_r\mathrm{<mo>=</mo>}\frac{n_{h_r}}{N_{h_r}}$ est la fraction de sondage dans la strate $h_r.$

En pratique, on utilise les estimations de la variance par linéarisation de Taylor et les estimations de la variance par la méthode du jackknife pour construire les intervalles de confiance $t$ de Wald avec $R-H$ degrés de liberté, où $R$ est le nombre d’UPE (ou le nombre de rééchantillonnages) et $H$ est le nombre de strates.

On montre facilement que l’estimateur de la variance par la méthode du jackknife équivaut algébriquement à l’estimateur linéarisé de Taylor pour les estimateurs linéaires du plan de sondage. En revanche, pour ce qui est des estimateurs non linéaires du plan de sondage, comme les coefficients de régression pour les modèles de régression à risques proportionnels, la méthode du jackknife tend à produire des estimations de la variance légèrement plus élevées que la méthode linéarisée de Taylor (Fuller, 2009).

Notons que si l’estimation de l’échantillon complet présente une vraisemblance monotone, il est très probable que la plupart des échantillons répétés présentent également des vraisemblances monotones. Il en résultera de nombreuses estimations par répliques « inutilisables ».

Les procédures d’analyse des données d’enquête dans SAS/STAT prennent en charge à la fois les méthodes d’estimation de la variance linéarisée de Taylor et d’estimation de la variance par répliques (Mukhopadhyay, An, Tobias et Watts, 2008).

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