Effet marginal ou conditionnel, quelle différence ?

Ce billet a pour objet de s’adresser à un public statistique moins aguerri que d’autres. Vous avez peut-être déjà entendu parler d’effet marginal ou conditionnel ou encore de modèle marginal ou conditionnel. Qu’est-ce que cela veut dire ou implique ?

Cette problématique s’applique aux régressions logistiques, aux régressions de Poisson, mais pas au modèle linéaire général. En bref, le problème se pose à partir du moment où la relation entre les variables explicatives et la variable à expliquer n’est pas linéaire. Cela s’illustre très bien avec un exemple. Considérons une population dans laquelle un risque (par exemple le décès à 30 jours) va toucher 52.5% de la population. Il existe un facteur pronostique majeur binaire, de prévalence 50%, divisant la population en deux sous-populations de taille égale : une sous-population à bas risque avec 10% de décès et une population à haut risque avec 95%. Comme les deux sous-populations représentent chacune 50% de la population globale, le risque général est bien de 0.50×0.10+0.50×0.95=52.5%.

Supposons qu’un traitement réduise le risque de 5% dans chacune des deux sous-populations. La sous-population à bas risque passe d’un risque de 10% à 5% et le haut risque de 95% passe à 90%. L’odds ratio protecteur du traitement, dans chacune des deux sous-populations est environ égal à 0.50. Pour être rigoureux, il est égal à 0.474 dans chacune des deux sous-populations. C’est l’odds ratio conditionnel au facteur pronostique.

Maintenant, on peut calculer que le risque général d’une population totalement exposée au traitement. La réduction de 5% du risque dans chacun des deux sous-groupes se traduit par une réduction de 5% du risque général qui passe de 52.5% à 47.5%. Cela peut aussi se calculer comme 0.50×0.05+0.50×0.90=47.5%. Ainsi, en calculant 52.5%-47.5% on retombe sur les 5% de réduction. La moyenne des différences (-5% dans chaque sous-population) est égale à la différence des moyennes (52.5% – 47.5%).

Grace à ce second calcul, nous pouvons calculer l’odds ratio protecteur du traitement sur la population générale. Le risque de 52.5% correspond à une cote de 0.525/(1-0.525)=1.105, qui est réduite à un risque de 47.5% soit une cote de 0.475/(1-0.475)=0.905. Au total, l’odds ratio est de 0.905/1.105 = 0.819. Cet odds ratio s’applique au pourcentage moyen de la population complète. C’est un odds ratio de la moyenne, ou odds ratio marginal.

Ainsi, nous avons noté que l’odds ratio de 0.819 est un odds ratio marginal alors que l’odds ratio de 0.475 est un odds ratio conditionnel au facteur pronostique, c’est-à-dire, s’appliquant à chacun des sous-groupes défini par ce facteur pronostique.

Que se passe-t-il lorsque les effets dans chacun des deux sous-groupes ne sont pas égaux ? Certains modèles, comme la régression logistique, reposent sur l’hypothèse de constance des effets conditionnels dans les sous-groupes. En cas d’écart à cette condition, en pratique, l’odds ratio calculé sera égal à une valeur intermédiaire, une sorte de moyenne de tous les odds ratio des sous-groupes, pondérée par la précision statistique de l’odds ratio dans ce sous-groupe. En bref, on peut à peu près obtenir cet odds ratio conditionnel « moyen » en calculant la moyenne des log-odds ratio pondérée par l’inverse de leur variance, puis en repassant à l’odds ratio par l’exponentielle. Selon la manière d’estimer le modèle, on obtiendra un résultat légèrement différent. Point important, il est à noter que les procédures de régression logistique implémentés par tous les logiciels fournissent toujours, en sortie, des odds ratio de chacune des covariable, conditionnels à l’ensemble des autres covariables.

On peut dire que l’odds ratio conditionnel s’applique séparément à chacune des observations, avec son risque de base propre (10% ou 95%), alors que l’odds ratio marginal est l’odds ratio qui s’applique à la moyenne de la population (risque de base 52.5%).

Quel est le problème des odds ratio conditionnels ? D’abord, le terme est incomplet… Conditionnel à quoi ? Un odds ratio conditionnel à l’âge, n’est pas le même qu’un odds ratio conditionnel au sexe, qui diffère encore de l’odds ratio conditionnel à l’âge et au sexe. Ensuite, plus le nombre de conditions est élevé et les facteurs pronostiques fortement liés à l’outcome, plus l’odds ratio conditionnel tend vers 0 ou l’infini. Dans l’exemple, l’odds ratio marginal était faible (0.819) alors que l’odds ratio conditionnel était nettement plus fort (0.475). Cela est explicable par le fait qu’une différence de 5% correspond à un odds ratio faible lorsqu’on est proche d’un risque de 50% (population générale) mais correspond à un odds ratio fort dès que le risque se rapproche de 0% ou 100%. Or, plus un modèle sera fortement prédictif, plus les sous-groupes auront des valeurs prédites extrêmes (proches de 0 et 100%) et plus les odds ratio conditionnels seront forts. Par ailleurs, mécaniquement, un modèle dans lequel on ajoute des covariables sera toujours plus fortement prédictif. Même si on ajoute des variables de corrélation nulle avec l’outcome, sur un échantillon fini, l’effet de cette variable sera légèrement positif ou négatif par hasard et cela améliorera la prédictivité (réduira la déviance) du modèle sur l’échantillon, avec pour effet de gonfler tous les odds ratio conditionnels.

Ainsi, deux odds ratio conditionnels à des choses différentes ne sont pas comparables. C’est pourquoi, au sens strict, on ne peut pas comparer un odds ratio dans un modèle avec une covariable à l’odds ratio de la même covariable dans un modèle avec une dizaine d’ajustements supplémentaires. De manière mécanique, l’odds ratio conditionnel à la dizaine de variables d’ajustement sera gonflé. En pratique, cela pose rarement problème en médecine, car les « facteurs pronostiques » sont généralement faiblement corrélés à l’outcome, conduisant à une inflation modeste des odds ratio conditionnels. Néanmoins, j’ai un beau cas d’école, pour lequel il existe une solution élégante. Le problème était d’évaluer les causes de la variance inter-établissement en terme de proportion de césariennes chez les femmes parturientes. La variance inter-établissement était exprimée comme l’odds ratio médian, interprétable comme la médiane des odds ratio entre deux établissements pris au hasard (le détail du calcul est un peu plus complexe et se basait sur la variance inter-établissement dans un modèle linéaire généralisé à effets mixtes, mais c’est approximativement ça). Or, les établissements n’ont pas la même population. Certains suivent des femmes enceintes à très bas risque alors que d’autres suivent des femmes à haut risque. La classification de Robson (doi:10.1016/j.jgyn.2015.02.001) définit 12 groupes de femmes, avec des risques extrêmement différents, allant de 2.1% pour les femmes multipares portant un singleton céphalique > 37 SA avec travail spontané, jusqu’à 82.4% pour les primipares avec présentation siège. Il était attendu que l’analyse brute montre un fort odds ratio inter-établissement médian, et qu’il s’atténue lors de l’ajustement sur le Robson. Au contraire, l’odds ratio inter-établissement médian a augmenté. Cela est explicable par le fait que des différences minimes dans les sous-groupes à très bas ou très haut risque (p.e. 2.1% -> 1% pour le sous-groupe à bas risque) s’expriment par des odds ratio conditionnels importants alors que l’effet marginal (20% -> 19%) s’exprime par un odds ratio bien plus modeste. Outre l’usage d’un modèle linéaire qui résout complètement le problème, on peut calculer dans chaque établissement un taux standardisé de césarienne en faisant la moyenne des taux des sous-groupes, pondérée par une fréquence de référence de ce sous-groupe de Robson. La fréquence de référence des sous-groupes peut se baser sur la littérature ou simplement être obtenu en poolant tous les établissements afin d’obtenir une population de référence représentative de l’ensemble des centres. Les odds ratio obtenus entre ces taux standardisés sont marginaux et sont alors comparables aux odds ratio d’un modèle non ajusté.

Enfin, il existe un problème à ne pas négliger : les modèles de régression logistiques dits « conditionnels » que l’on utilise en cas d’appariement, comme, par exemple dans les études en cross-over fournissent un odds ratio ininterprétable. Pour faire bref, ils fournissent une probabilité conditionnelle à une condition inobservable. Cela s’illustre bien en considérant le cas d’une maladie chronique stable pour laquelle on essaye un traitement symptomatique (p.e. traitement de la douleur dans la polyarthrite rhumatoïde) et pour laquelle un essai en cross-over est pertinent. Considérons pour critère de jugement, la réponse binaire au traitement, défini par une amélioration de la douleur. Considérons l’usage de deux traitements de nature très proche et tous deux assez inefficaces. On peut imaginer que 90% des sujets ne répondront à aucun des deux traitements, 8% des sujets répondront aux deux traitements, 1.5% des sujets répondront au traitement innovant mais pas au traitement de référence et 0.5% des sujets répondront au traitement de référence mais pas au traitement innovant. Dans ce contexte, les taux de réponses avec le traitement de référence et le traitement innovant sont respectivement de 8.5% et 9.5%, soit un différence absolue de pourcentage de 1% (nombre de sujets à traiter pour obtenir une réponse supplémentaire = 100) ou un odds ratio marginal à 1.13. L’odds ratio conditionnel au patient est égal à 3.0. En effet, cet odds ratio conditionnel au patient est calculé seulement sur les 2% de paires de mesures discordantes. Les 80% de paires concordantes négatives ne participent pas à l’estimation de l’odds ratio car les chances de réponses sont de 0% dans cette paire et tous les odds ratio ont la même vraisemblance . De même pour 8% de patients répondant aux deux traitements. La statistique est seulement calculée comme le rapport entre les discordances favorables au traitement innovant (1.5%) et les discordances favorables au traitement de référence (0.5%). Le rapport 1.5%/0.5% est égal à 3.0, l’odds ratio conditionnel au patient, ou encore, pour être plus clair, conditionnel au fait que le patient ait une réponse différente aux deux traitements. Cet odds ratio à 3.0 est extrêmement élevé par rapport à l’odds ratio marginal car la corrélation intra-patient est forte, ou, en d’autres termes, le facteur patient est fortement pronostique ! Le problème, c’est que ce facteur patient est inobservable tant qu’on n’a pas essayé les deux traitements. Si on considère que les patients se divisent en trois catégories : les rouges (80%), qui ne répondent à aucun des deux traitements, les verts (8%) qui répondent aux deux traitements et les oranges (2%) qui répondent à l’un mais pas à l’autre, alors l’odds ratio fourni s’applique au sous-groupe des patients oranges. Il est conditionnel à la couleur orange. Si la couleur était une variable clinique identifiable, alors il serait aisé d’orienter la prescription. Les patients rouges n’auraient aucun de ces deux traitements, les patients verts auraient l’un ou l’autre, à la préférence du médecin ou du patient et les patients oranges auraient en priorité le traitement innovant avec trois fois plus de chances de réponse que si on leur donnait le traitement de référence. Problème majeur : la couleur est inobservable, sauf à réaliser une période d’essai cross-over pour le patient qui conduirait à identifier d’emblée le sens de la différence et donc, à l’inutilité de la connaissance de cet odds ratio. En bref, l’odds ratio est conditionnel à une donnée inconnue.

Le problème décrit avec la régression logistique conditionnelle s’applique aussi aux régressions logistiques à effets mixtes dont les effets sont conditionnels aux effets aléatoires qui reposent sur des variables inobservables. Ces considérations n’existent pas dans les modèles linéaires à effets fixes ou à effets mixtes gaussiens parce que les effets conditionnels dans ces modèles sont égaux aux effets marginaux la fonction de lien étant l’identité. C’est pourquoi l’essai en cross-over imaginaire cité ci-dessus pourrait être analysé avec une estimation de la différence absolue de chances de réponse par la méthode de Student sur séries appariées. Cette méthode est asymptotiquement correcte et fournit un résultat pertinent et indépendant du degré de corrélation intra-paire.

Au final, rien n’est plus simple et interprétable que les résultats de modèles linéaires gaussiens. Les odds ratio conditionnels risquent d’être incorrectement interprétés comme des odds ratio marginaux alors qu’ils sont toujours gonflés de manière plus ou moins extrême selon que la condition explique plus ou moins fortement la variance de l’outcome. Les modèles logistiques dits ‘conditionnels’ sont à éviter car fournissent des résultats théoriquement ininterprétables, en pratique, fournissent des odds ratio interprétés à tort comme marginaux, et donc plus ou moins fortement biaisés.

Modèle multivarié et essai clinique randomisé

Résumé

Ce billet discute de l’ajustement dans les essais cliniques en randomisation individuelle. Il montre que l’ajustement statistique fait peu gagner en puissance, permet parfois de rectifier un défaut de randomisation, mais pose aussi des soucis non négligeables :

  1. P-hacking plus ou moins important
  2. Erreurs de manipulation du logiciel
  3. Résultats théoriquement non biaisés mais ininterprétables (odds ratio conditionnel à une variable inobservée), en pratique interprétables mais biaisés (car interprétés comme odds ratio marginal ou conditionnel aux variables observées)
  4. Résultats mal présentés (p.e. odds ratio conditionnel plutôt que différence absolue de risques)
  5. Les hypothèses faites par les modèles (absences d’interactions, linéarité) sont fausses, et cela peut conduire à un biais dans les estimations

Une partie de ces problèmes disparaît ou s’atténue lorsqu’on utilise exclusivement le modèle linéaire gaussien, qu’on ajuste seulement sur la covariable pronostique majeure, qu’on s’assure que le protocole soit publié avec un paragraphe statistique très exhaustif.

Je conseille d’avoir conscience de ces problèmes afin d’utiliser avec parcimonie, et en toute connaissance de cause, les ajustements dans les essais cliniques en randomisation individuelle.

Introduction

Dans un essai clinique randomisé bien mené, il n’y a pas de facteur de confusion car l’allocation des traitements est aléatoire et donc indépendante de tout facteur pronostique ou pas. Un ajustement statistique sur des facteurs pronostiques est possible, avec ou sans stratification de la randomisation. Cet ajustement réduit les fluctuations d’échantillonnage explicables par les déséquilibres aléatoires entre les groupes sur ces facteurs pronostiques, ou, dans le cadre de la randomisation stratifiée, évite une surestimation des fluctuations d’échantillonnages, en excluant du calcul de la variance résiduelle, la variance ou plutôt la non-variance des facteurs pronostiques. Au final, cela augmente la puissance sans modifier le risque alpha. Comme dans l’évaluation de toute méthode statistique, deux questions doivent être posées : que gagne-t-on ? que perd-on ?

Que gagne-t-on ?

Le gain de puissance est maximal lorsque le facteur pronostique expliquant une grande partie de la variance de l’outcome. En pratique, cela veut dire que le facteur pronostic est fortement lié à l’outcome et a une grande variance. Pour un facteur pronostique binaire, l’idéal est une prévalence de 50% (forte variance) et que selon sa présence ou non, l’outcome soit très favorable ou défavorable.

Nous prendrons comme référence, un cas d’école dans lequel le facteur pronostique est majeur et a une forte variance : le stade selon l’American Joint Committee on Cancer (AJCC) du cancer du colon, communément appelé stade TNM. Nous prenons comme référence la version 6 (O’Connell J, Maggard M, Ko C. Colon Cancer Survival Rates With the New American Joint Committee on Cancer Sixth Edition Staging. Journal of National Cancer Institute, Vol. 96, No. 19, October 6, 2004) car elle est en accès libre et dispose de tables de mortalité et de prévalences précises contrairement aux versions 7 et 8. La figure ci-dessous montre la distribution jointe du stade et de la mortalité à 5 ans.

Un facteur pronostique capable de séparer 8.1% de mortalité de 93.2%, ainsi que de nombreux risques intermédiaires, avec une variance aussi grande n’est pas commun dans les essais cliniques. Il faut noter qu’un essai clinique n’inclura pas des stades I et des stades IV. Le pronostic et le traitement de référence dépendant beaucoup du stade, un essai clinique inclura un seul stade (exemples : https://dx.doi.org/10.1056/NEJMoa1713709, https://dx.doi.org/10.1001/jamaoncol.2019.6486). Autrement, mélanger des patients ayant des stades très différents ouvre la porte à des interactions et la sélection d’un traitement bénéfique à un sous-groupe (p.e. stade III) mais nocif à un autre (p.e. stade II), sans que la distinction ne soit faite entre les deux stades. En bref, un tel facteur pronostic majeur ne devrait pas exister à l’intérieur d’un essai clinique. Considérons dans un premier scenario qu’il soit possible d’inclure des stades I (tumeur localisée ne dépassant pas la muqueuse) et des stades IV (tumeur avec métastases à distance) dans un même essai clinique.

Nous ne nous intéresserons pas aux modèles de survie basés sur les risques proportionnels, qui ont leur propres problèmes, mais plutôt considèrerons la mortalité binaire à 5 ans.

Le 1er scenario est basé sur une puissance à 80%, un risque alpha unilatéral à 2.5% (ou bilatéral à 5%) dans un modèle ajusté sur le stade tumoral (7 modalités -> 6 variables binaires recodées) en randomisation simple 1:1 dans un modèle linéaire gaussien pour une intervention correspondant à un odds ratio conditionnel à 0.50 (facteur protecteur) sur la mortalité à 5 ans. L’omission d’un ajustement est à l’origine d’une inflation de l’erreur type (standard error) par un facteur 1.32, correspondant à une réduction de la puissance de 80% à 56.4%. Pour atteindre la même puissance à 80%, il faut multiplier par 1.75 le nombre de sujets, passant de 519 sujets à 907 sujets. Le gain en terme de puissance est peu dépendant de l’effet du traitement. Un odds ratio à 0.25 conduit à une perte de puissance de 80% à 57.1% lors du non-ajustement alors qu’un odds ratio à 0.95 conduit à une perte de puissance de 80% à 57.4% lors du non-ajustement. Cela est dû au fait que la perte de puissance dépend seulement de la variance expliquée par le facteur pronostique, qui dépend de la prévalence et de la moyenne des outcomes les deux groupes étant poolés. Cette moyenne des outcomes dépend un peu de l’effet du traitement.

Le 1er scenario montre un intérêt important de l’ajustement mais n’est pas réaliste du tout. Aucun essai clinique n’inclurait des stades I à IV. Le 2ème scenario, encore loin de la réalité, est basé sur l’inclusion des stades I à III (tumeur très localisée, jusqu’à tumeur avec envahissement local important touchant les organes voisins et métastases ganglionnaires régionales nombreuses). Dans ce cas là, le non-ajustement (avec odds ratio à 0.50) conduit à une inflation de l’erreur type par 1.061, une réduction de la puissance de 80% à 75.2% ou une augmentation du nombre de sujets nécessaires de 12.5%. Comme dans le 1er scenario, l’effet du traitement n’a pas vraiment d’importance.

Le 3ème scenario est identique au second même une régression logistique est employée plutôt qu’une régression linéaire. La perte de puissance est presque la même que dans le modèle linéaire, à une différence absolue de puissance de 0.98% près. Par contre, le modèle logistique risque d’être mal interprété. En effet, le bénéfice de traiter tous les sujets avec le traitement innovant correspond à une réduction de la mortalité à 5 ans de 38.4% à 29.6% correspondant à une différence absolue de risque de 8.8% (meilleur manière d’en quantifier le bénéfice) ou à un odds ratio marginal de 0.675. L’odds ratio conditionnel au stade, par contre, est égal à 0.50 (soit 0.74 fois l’odds ratio marginal). C’est cet odds ratio conditionnel que les logiciels statistiques fournissent en sortie. Plus l’ajustement porte sur un facteur pronostique pertinent, plus l’odds ratio conditionnel s’écarte de l’odds ratio marginal. En l’absence d’ajustement, la régression logistique fournit l’odds ratio marginal, plus simple à interpréter car indépendant du nombre et de la nature des variables d’ajustement.

Le 4ème scenario, plus réaliste, consiste en l’ajustement sur le stade a, b ou c chez des patients tous en stade III, en conservant les autres paramètres (odds ratio conditionnel à 0.50, puissance à 80%, risque alpha à 2.5% unilatéral), dans un modèle linéaire gaussien. L’inflation de la variance due au non-ajustement est de 1.025 avec une perte de puissance de 1.97% (80% -> 78.03) ou un nombre de sujets nécessaires augmenté de 5.1%.

En bref, l’ajustement sur un facteur pronostique important et de forte variance, dans des conditions réalistes, peut faire gagner 2% de puissance soit 5% de sujets à recruter en moins. Il paraît déjà impossible de gagner 5% de puissance car il faudrait une hétérogénéité majeure de la population (stades I à III mélangés), ce qui conduirait à une réduction du nombre de sujets nécessaires d’environ 11%. On peut donc conclure que le gain, en nombre de sujets nécessaires sera toujours inférieur à 10%, et sera souvent situé en dessous de 5%. La présence de multiples facteurs pronostiques fournit rarement une capacité pronostique supérieure à celle du stade dans les cancers du colon. Par exemple, le Karnofsky n’est pas pris en compte dans le stade AJCC, probablement parce qu’il n’influence pas de manière majeure le pronostic.

Un second bénéfice, non négligeable de l’ajustement, c’est l’adaptation aux situations où l’essai clinique randomisé n’a pas été bien mené telle qu’une mauvaise randomisation ce qui peut arriver, par exemple, dans un essai en ouvert avec une procédure de randomisation par blocs de permutation de taille 4 stratifiée sur le centre, permettant rapidement au clinicien de deviner dans quel groupe un patient sera randomisé. Cela peut aussi parfois servir s’il y a eu un biais d’attrition différentiel. L’ajustement peut, dans une certaine mesure, rectifier le biais.

Que perd-on ?

L’ajustement dans les essais cliniques randomisés n’est pas dénué de risques. Même si les essais cliniques sont très souvent enregistrés dans des bases de données comme ClinicalTrials ou EudraCT, le résumé de protocole disponible précise le critère de jugement mais rarement l’analyse statistique. Les protocoles complets détaillant l’analyse statistique ne sont pas toujours disponibles. Si le protocole complet n’est pas disponible les ajustements peuvent être un facteur conscient ou inconscient de P-hacking. Il est fréquent de lancer un très grand nombre d’analyses de sensibilité dans les rapports d’analyse statistique et l’analyse qui permettra le passage du petit p en dessous de 0.05 risque d’être choisie comme analyse principale.

Même un protocole détaillant la méthode, laisse généralement des degrés de liberté.

Exemple : l’analyse principale sera faite en intention de traiter modifiée (exclusion des sujets randomisés sortis d’étude ou décédés avant d’avoir reçu la 1ère cure de traitement) et sera réalisée dans un modèle de Cox expliquant la survie sans progression ajustée sur le centre, l’âge, le Karnofsky, le stade tumoral et le type histologique. Les sujets perdus de vue seront censurés à la date des dernières nouvelles. Le hazard ratio de l’effet traitement sera comparé à la valeur 1 par un test bilatéral au seuil de significativité 5%.

Cela semble clair et précis, mais il manque un grand nombre de détails :

  1. Par quelle méthode les ex-aequo seront-ils gérés : Efron, Breslow ou méthode exacte ?
  2. Certains ajustements sont-ils réalisés par stratification (modèle de Cox conditionnel) ? Cela est parfois fait pour l’effet centre.
  3. Le petit p est-il calculé par le test du rapport de vraisemblance, celui du score ou celui de Wald, ou encore un autre estimateur ?
  4. Comment seront gérées les données manquantes sur les covariables, notamment le Karnofsky ?
  5. Comment seront recodées les variables. L’âge, par exemple, peut être introduit sous forme d’une unique variable quantitative (effet log-linéaire), sous forme de polynôme, ou découpé en catégories à seuils prédéfinis (p.e. 18-60, 60-75, >75 ans) ou découpé en quantiles (quintiles, quartiles, tertiles).
  6. Quel est le niveau de détail du stade ? Distingue-t-on le stade IIa du stade IIb ou regroupe-t-on ces deux modalités en stade II ?
  7. Quel est le niveau de détail du type histologique ? En effet, on peut toujours sous-typer très précisément avec la biologie moléculaire, ou pas…
  8. La variable de stade fait-elle référence au cTNM, au pTNM ? Est-ce le stade pré-traitement ou post-traitement néo-adjuvant ?
  9. Comment sont gérées les valeurs aberrantes dans les covariables découvertes après le lever d’aveugle ?
  10. Quelles conditions de validité seront testées ? Y aura-t-il des tests d’interaction ? Si oui, que fera-t-on si certains sont significatifs ?

Pour l’item 10, la réponse évidente (de mon point de vue), est la non prise en compte des interactions, mais il n’est pas sûr que ce soit le comportement de tous les statisticiens. De même, pour l’item 9, la réponse évidente est l’imputation de la valeur aberrante par la même méthode de gestion des données manquantes que les autres, évitant ainsi un biais de classement différentiel, mais j’ai souvent observé dans ma pratique de biostatisticien, le comportement correspondant à « corriger » la donnée à partir du dossier médical, en gardant la trace de la modification dans le script d’analyse.

En analyse non ajustée, il existe déjà un espace de liberté permettant du P-hacking (p.e. plusieurs estimateurs existent, la gestion des ex-aequo se pose toujours), mais ces possibilités explosent complètement en cas d’analyse ajustée. Enfin, il est toujours possible de ne pas respecter du tout le protocole, sans que ce soit pour autant visible dans l’article. J’ai déjà observé un essai clinique randomisé dont les résultats sont publiés dans le Lancet pour lequel le protocole spécifiait que l’analyse principale serait faite dans une régression logistique ajustée sur la sévérité de la maladie à baseline et d’autres variables qui ressortiraient dans l’analyse univariée. Au final, plusieurs techniques de pas-à-pas ont été utilisées et l’article a été publié avec un modèle log-binomial non ajusté sur la sévérité de la maladie à baseline mais ajusté sur un score clinique à baseline qui était ressorti dans les analyses pas-à-pas (basées sur la P-valeur puis sur l’AIC). Cet exemple ne prouve pas que ces écarts au protocole sont fréquents, mais qu’ils sont possibles. Comme le protocole n’a pas été publié, cela est invisible.

Au delà du P-hacking conscient ou inconscient, il existe le risque d’erreur d’analyse. La manipulation des outils statistiques n’est pas évidente et il est toujours possible de se tromper dans la programmation d’une imputation, d’un recodage de variable ou de l’estimation d’un modèle. Plus le nombre de variables impliquées est grand et plus le risque d’erreur est, a priori, élevé. Je ne dispose pas de données à ce jour pour évaluer la fréquence ou l’ampleur de ce problème.

Il apparaît aussi le problème du choix du modèle, pour lequel il est aisé de faire le mauvais choix. Sur un critère de jugement binaire, appliquera-t-on un modèle de régression logistique (logit-binomial), un modèle log-binomial, un modèle identité-binomial ou un modèle linéaire gaussien ? Si on ajuste sur l’effet centre, va-t-on utiliser les équations d’estimation généralisées (GEE), un modèle linéaire généralisé à effets mixtes (GLMM) ou un modèle linéaire généralisé à effet fixe ? Le problème, c’est que l’analyse univariée fournira aisément un résultat pertinent, alors que peu des modèles sus-cités fournissent des résultats aussi robustes. Une régression logistique fournit un odds ratio conditionnel ininterprétable dans un GLMM car conditionnel à des variables inobservables. Une régression logistique ou un modèle log-binomial fournit un odds ratio conditionnel rarement converti en effet marginal. Comme vu au-dessus, il est « gonflé » par rapport à l’effet marginal. Ensuite, la régression logistique comme la log-binomiale fournissent des statistiques relatives (odds ratio ou rapports de proportions) qui reflètent bien mal le bénéfice réel du traitement dans un essai clinique randomisé. En effet, un risque relatif à 0.80 sur un risque touchant la moitié des sujets correspond à une réduction de 10% du risque absolu, soit seulement 10 sujets à traiter pour éviter un événement. Par contre, le même risque relatif appliqué à un risque extrêmement rare tel que 0.1%, obligera à traiter 5000 personnes pour éviter un événement. C’est pour cela, par exemple, qu’une antibioproxphylaxie dans des chirurgies à très haut risque infectieux se justifie alors qu’elle serait futile dans des opérations à très bas risque.

Le modèle identité-binomial semble pertinent, évaluant directement une différence absolue de risque, mais celui-ci montre ses limites lorsque des interactions apparaissent dans le modèle. Comment estimer le bénéfice net d’un traitement lorsque la différence absolue de risque diffère selon le sous-groupe ? En présence d’une interaction quantitative, c’est l’effet marginal qui nous permet d’évaluer le bénéfice net. Si 10% des sujets (mauvais pronostic) voient leur risque d’évolution défavorable baisser de 30% alors que 90% des sujets (bon pronostic) voient leur risque d’évolution défavorable baisser de 20% alors en donnant le traitement innovant à l’ensemble des sujets ont réduit de 0.10×0.30 + 0.90×0.20 = 21%. Un modèle linéaire non ajusté estimera, sans biais, les 21% de bénéfice. Un modèle linéaire gaussien ajusté estimera aussi, sans biais, les 21% de bénéfice car l’homoscédasticité supposée par le modèle conduira à une statistique ajustée égale à la des sous-groupes pondérée par l’effectif du sous-groupe. Cela est une propriété de la méthode d’estimation des moindres carrés. Le modèle identité-binomial, par contre, supposera que la variance est plus faible pour les pourcentages proches de 0% ou de 100% et donnera donc un poids plus fort au sous-groupe correspondant. En caricaturant volontairement, si la population est composée d’un groupe de 80% de patients ayant un risque à 99% (réduit de à 89% par le traitement, soit -10%) et 20% de patients ayant un risque à 50% (réduit à 49% par le traitement, soit -1%), alors le modèle identité-binomial ajusté sur le facteur pronostique estimera le bénéfice du traitement à -6.2% (espérance de l’estimateur) contre -2.8% pour un modèle identité-binomial non ajusté ou un modèle linéaire gaussien ajusté ou non ajusté. Le vrai bénéfice, si on donne le traitement à tout le monde, correspond à cette seconde estimation.

Il est illusoire de penser que les interactions quantitatives n’existent pas. Déjà, par principe, s’il n’existe pas d’interaction quantitative dans un modèle logistique, alors il en existera dans le modèle linéaire et vice versa. Les deux modèles sont fondamentalement contradictoires (sauf s’ils contiennent tous deux tous les termes d’interactions possibles). Les autres modèles sont aussi contradictoires : log-binomial, probit, etc. La stratégie consistant à tester les interactions et à accepter l’hypothèse nulle d’absence d’interaction par défaut de preuve du contraire, serait inepte à un bayésien qui sait que la probabilité a priori d’interaction est à peu près égale à 100%. C’est pourquoi, il est préférable de construire des statistiques « robustes » aux interactions quantitatives. Pour cela, il est possible d’utiliser un modèle linéaire gaussien ou identité-binomial non ajusté ou un modèle linéaire gaussien ajusté. Il est aussi possible de se baser sur un modèle AVEC interactions dès le départ, puis de reconstruire la statistique marginale en calculant la moyenne des effets dans chacun des sous-groupe pondérée par la taille du sous-groupe. Comme la covariance entre les sous-groupes devient nulle lorsqu’on utilise un modèle avec des interactions complètes, le calcul n’est pas super compliqué. Le choix du modèle n’a alors plus d’importance puisqu’on estime directement les effets moyens en sous-groupes. En présence de variable pronostique quantitative, il n’y a pas de solution de codage parfaite parce qu’on ne peut pas modéliser la nature de la relation de manière exacte. Le problème disparaît si on n’ajuste pas du tout ou si on découpe la variable quantitative en tranches (ajustement partiel).

En cas d’interaction qualitative, les problèmes décrits ci-dessus sont exacerbés. C’est probablement un phénomène rare, heureusement, mais il est préférable de se comporter aussi bien que possible s’il advenait. Dans le cas d’une interaction qualitative, le modèle identité-binomial est susceptible de montrer un effet opposé à l’effet marginal réel ! Les interactions qualitatives sont problématiques dans tous les cas parce que le traitement est bénéfique à un sous-groupe identifiable de patients alors qu’il est nocif à un autre sous-groupe identifiable de patients. Malheureusement, la puissance statistique nécessaire à la détection de ce phénomène est loin d’être satisfaisante dans la plupart des études. Cela est d’autant plus vrai que les sous-groupes sont déséquilibrés. C’est pourquoi on évitera, dès le départ, par des critères d’inclusion judicieux, d’avoir une population trop hétérogène avec une probabilité d’interaction forte. Cela justifie, par exemple, qu’on évite de mélanger des enfants et des adultes dans les évaluations de traitements psychiatriques. Néanmoins, dans le cas où on se fait avoir, c’est-à-dire, où il existe une interaction qualitative sur une variable qu’on n’avait pas anticipé et pour laquelle la puissance est totalement insuffisante pour retrouver une tendance dans une analyse en sous-groupes, alors on devra « vivre » avec et donner le traitement innovant à tout le monde. C’est la même chose avec les variables d’interaction inconnues (p.e. dosage biologique encore inconnu) : en l’absence de connaissance, on évalue le bénéfice global. C’est là que l’effet marginal directement estimé par le modèle linéaire gaussien reflètera le bénéfice pragmatique. Empiriquement, il semblerait que le modèle de régression logistique ne puisse pas conduire à l’erreur de 3ème espèce (effet significatif dans le sens opposé à l’effet réel) là où le probit-binomial, l’identité-binomial et le log-binomial peuvent. Il semblerait qu’en présence d’interaction, même si elle n’est pas prise en compte dans le modèle logistique, l’effet marginal prédit par le modèle logistique est correct (identique au modèle linéaire ou au modèle non ajusté) dans le cas d’un ajustement sur une covariable catégorielle.

Il existe aussi le problème de non convergence très fréquent dans le modèle log-binomial lorsqu’on ajuste sur une variable quantitative ou sur de nombreuses variables binaires. Cela est dû au fait que le modèle est grossièrement faux, conduisant à des prédictions dépassant le risque de 1. Il y a de nombreuses interactions quantitatives ou défauts de log-linéarité dans ces modèles et comme vu au-dessus, il n’est plus possible d’estimer un effet conditionnel ou marginal de manière fiable. Par ailleurs, que faire en cas de non convergence ? Fournir un effet non ajusté ? Cela paraît raisonnable, mais la procédure en deux étapes entraîne une erreur statistique mal caractérisée. C’est par exemple, lorsque l’effet d’une covariable est surestimé, que le défaut de convergence apparaît, et dans ce cas précis, on supprime l’ajustement, supposant forçant son effet est à zéro et créant une fluctuation d’échantillonnage discrète sur la statistique principale. En bref, en ajoutant un seul sujet, on peut changer l’estimation principale de manière principale. Sous l’hypothèse d’absence d’effet de la variable principale, cela entraîne une erreur mais pas un biais car ça n’a aucune raison d’avantager un groupe que l’autre. Sous l’hypothèse d’existence d’un effet, cela pourrait peut-être engendrer un biais minime, car il y a à peu près une chance sur deux que cela fasse fluctuer la statistique principale dans un sens et une chance sur deux que ça la fasse fluctuer dans l’autre étant donné que la corrélation entre la covariable et la variable principale est totalement aléatoire de moyenne que la variable

Il existe aussi le problème de non convergence très fréquent dans le modèle log-binomial lorsqu’on ajuste sur une variable quantitative ou sur de nombreuses variables binaires. Cela est dû au fait que le modèle est grossièrement faux, conduisant à des prédictions dépassant le risque de 1. Il y a de nombreuses interactions quantitatives ou défauts de log-linéarité dans ces modèles et comme vu au-dessus, il n’est plus possible d’estimer un effet conditionnel ou marginal de manière fiable. Par ailleurs, que faire en cas de non convergence ? Fournir un effet non ajusté ? Cela paraît raisonnable, mais la procédure en deux étapes entraîne une erreur statistique mal caractérisée. C’est par exemple, lorsque l’effet d’une covariable est surestimé, que le défaut de convergence apparaît, et dans ce cas précis, on supprime l’ajustement, supposant forçant son effet est à zéro et créant une fluctuation d’échantillonnage discrète sur la statistique principale. En bref, en ajoutant un seul sujet, on peut changer l’estimation principale de manière principale. Sous l’hypothèse d’absence d’effet de la variable principale, cela entraîne une erreur mais pas un biais car ça n’a aucune raison d’avantager un groupe que l’autre. Sous l’hypothèse d’existence d’un effet, cela pourrait peut-être engendrer un biais minime, car il y a à peu près une chance sur deux que cela fasse fluctuer la statistique principale dans un sens et une chance sur deux que ça la fasse fluctuer dans l’autre étant donné que la corrélation entre la covariable et la variable principale est totalement aléatoire et d’espérance nulle. Difficile de dire si l’erreur aléatoire augmente un peu le risque alpha ou le diminue un peu. En tout cas, l’erreur aléatoire ne suit pas une loi normale.

Une autre option, en cas de non-convergence d’un modèle log-binomial, c’est d’utiliser un modèle tout autre à ce moment là, tel qu’un modèle logistique. Cela est problématique parce que la statistique fournie n’est plus comparable. Les fluctuations d’échantillonnage fréquentistes deviennent presque ininterprétables puisque la base de la statistique fréquentiste, c’est qu’un même protocole conduit à une même statistique.

Le modèle identité-binomial a des défauts de convergence aussi, même s’ils sont beaucoup plus rares que dans le modèle log-binomial. L’usage d’un modèle linéaire gaussien avec des ajustements parcimonieux (une ou deux variables pronostiques majeures, les autres étant généralement futiles) évite le problème.

Le problème de non-convergence est aussi une problématique à prendre en compte dans le protocole, qui doit préciser l’usage de modèles de secours. Autrement, cela pourrait être un facteur de P-hacking ou de non reproductibilité des analyses.

Synthèse

Nous avons vu qu’un ajustement statistique dans un essai clinique randomisé permet de recruter 5% de patients en moins en cas de facteur pronostic important et que ce chiffre peut monter jusqu’à 10% dans les conditions ou un facteur pronostic absolument majeur existe (en mélangeant les stades I à III du cancer du colon dans un même essai clinique par exemple). Il existe aussi un bénéfice de rectification des biais dû à une étude mal conduite, telle qu’une triche avec la randomisation.

Par contre, l’ajustement pose des problèmes de P-hacking et a de nombreuses limites d’interprétation à moins qu’il ne soit réalisé de manière très rigoureuse:

  1. Précision exacte de la méthode de traitement des données manquantes des covariables
  2. Précision du paramétrage exact du modèle, tel que l’estimateur
  3. Précision du codage exact des covariables rentrées dans le modèle
  4. Usage d’un modèle linéaire gaussien pour une variable binaire plutôt que les régressions logistiques, log-binomiales ou identité-binomiales
  5. Publication du protocole avec le plan d’analyse statistique le plus tôt possible !

Cela est plus facile à réaliser si vous ajustez sur peu de covariables et vous concentrez seulement sur une ou deux covariables pronostiques majeures, pour laquelle la qualité de la donnée est bonne (donnée de routine bien renseignée). S’il y a plusieurs covariables, ajoutez les termes d’interaction afin d’éviter les problèmes de fausses hypothèses.

Si vous êtes méthodologiste, je vous invite à être extrêmement explicite dans le protocole même pour un modèle non ajusté. Je vous conseille de ne faire un ajustement qu’à condition (i) de vous assurer que le protocole sera publié et (ii) qu’il existe un facteur pronostique majeur.

Si vous lisez un article, je vous conseille de toujours essayer d’obtenir le protocole, et de vérifier, s’il existe, qu’il a été respecté, notamment sur le choix du critère de jugement principal et le choix de l’analyse principale. Si le protocole n’est pas publié ou qu’il est incomplet, méfiez vous des analyses ajustées.

Randomisation simple ou pas

Les randomisations par blocs de permutation, de tailles fixe ou aléatoire, sont très utilisées dans les essais cliniques randomisés, qu’ils soient en aveugles ou en ouverts. Il arrive souvent que la randomisation soit stratifiée sur un ou plusieurs facteurs, tels que le centre ou des facteurs pronostiques, multipliant ainsi les listes de randomisation. D’un point de vue théorique, ces méthodes réduisent les déséquilibres en effectif entre les groupes (et sous-groupes) et améliorent la stabilité des estimateurs, conduisant à une meilleure puissance statistique. Malheureusement, il semblerait que le gain de puissance soit minime et qu’il y ait un risque de biais d’allocation, dans les études en ouvert.

Nous décrirons d’abord, un point de vue épistémologique, puis fournirons les détails statistiques permettant d’apprécier le gain de puissance et la perte de validité

Point de vue épistémologique

Avant d’utiliser la randomisation, d’autres moyens d’allocation pseudo-aléatoire de traitements avaient été utilisés, notamment l’inclusion dans le groupe A les jours pairs du calendrier et l’inclusion dans le groupe B les jours impairs. Cette méthode, semble, à première vue, fournir des groupes comparables, puisqu’il n’y a aucune raison que les caractéristiques du patient puissent être corrélées à la parité du jour. Néanmoins, cette méthode avait le défaut d’être prévisible. En situation où de grands espoirs sont portés sur le nouveau traitement, notamment lorsque le traitement de référence est peu efficace, la maladie sévère, et les bénéfices potentiellement perçus du traitement innovants sont importants, alors l’investigateur est susceptible d’utiliser cette information pour forcer le pseudo-hasard dans le sens qu’il souhaite. Souhaitant faire bénéficier du traitement innovant le patient dont le pronostic sous traitement de référence est péjoratif, il pourra décaler la consultation afin de s’assurer que celle-ci tombe un jour favorable. Dès lors qu’une information partielle ou totale sur le groupe dans lequel le patient se retrouverait existe avant l’inclusion, le destin peut être modifié par la volonté humaine. Dans le pire des cas, l’inclusion est différée et une nouvelle consultation est programmée. Si ce destin est modifié de manière répétée, un déséquilibre apparaît entre les deux groupes : il ne sont plus comparables. Le biais d’indication des études observationnelles ré-apparaît partiellement ou pleinement.

Le seul moyen d’éviter ce biais est la parfaite imprévisibilité de la séquence d’allocation. La randomisation, c’est-à-dire l’allocation au hasard, garantit cette imprévisibilité car le hasard n’est corrélé à aucune variable observée ou inobservée. La randomisation par blocs, ne répond pas à cette définition de hasard. Elle a une auto-corrélation négative et est donc prévisible dès lors que la randomisation est ouverte. Un investigateur incluant 4 patients d’affilée dans le même groupe peu parier que le prochain patient sera alloué dans l’autre groupe. L’algorithme est tellement simple, que son exploitation peut être inconsciente.

Gain de puissance

Le code R ci-dessous calcule la perte de puissance due à un déséquilibre des groupes en randomisation 1:1 avec une unique liste de randomisation:

power_for_n1=function(power_for_n, sig.level, n1, N) {
# On recrute N sujets en randomisation simple
# On observe n1 sujets dans le groupe A
# On veut connaître la puissance (conditionnelle à n1)
df = (N-2)
# espérance de la statistique T de Student pour n = N/2
# (randomisation par blocs -> équilibre parfait des groupes)
T_for_n = qt(power_for_n,df=df) + qt(1-sig.level/2,df=df)

# taille de n2 sachant que n1+n2 = N
n2 = N-n1
n  = N/2
# espérance de la statistique T de Student
# pour la randomisation simple (conditionnelle à n1)
T_for_n1 = T_for_n * sqrt((1/n + 1/n)/(1/n1 + 1/n2))

# puissance pour la randomisation simple (conditionnelle à n1)
pt(T_for_n1 - qt(1-sig.level/2,df=df), df=df)
}
# n1 suit une loi binomiale
power_for_nrandom=function(power_for_n, sig.level=0.05, N=100) {
x = 1:(N-1) # valeurs possibles pour n1
d = dbinom(x, N, 0.50) # probabilités associées
p = power_for_n1(power_for_n, sig.level, x, N) # puissances conditionnelles
sum (d * p) # puissance inconditionnelle
}

On peut ainsi calculer qu’une randomisation simple 1:1 avec une seule liste, fait passer une puissance de 80% à 79.8% (réduction de 0.2%) pour une étude avec 200 sujets en tout (~100 par groupe).

Biais d’allocation

Un coût non négligeable existe dans les études en ouvert : le groupe auquel va appartenir le prochain patient potentiellement incluable est partiellement prévisible.
Par exemple, dans une randomisation avec blocs de permutation de taille 4, il ne peut jamais y avoir plus de 4 patients d’affilée dans le même groupe. Cela est possible qu’un bloc est ‘0011’ avec le bloc suivant à ‘1100’. Dans une étude monocentrique, un investigateur qui a inclus 4 fois d’affilée des patients dans le même groupe sait alors dans quel groupe le prochain patient sera affecté. Dans une moindre mesure, si trois patients sont affectés d’affilée dans le même groupe, il y a 70.8% de chances que le prochain patient soit affecté à l’autre groupe. Le code ci-dessous simule ce scenario:

set.seed(2020)
v=as.vector(replicate(1000000,sample(c(1,1,0,0),replace=F)))

count_succ=function(v) {
	w = rep(0, length(v))
	for(i in 2:length(v)) {
		if (v[i] == 1) {
			w[i] = w[i-1]+1
		}
	}
	w
}


w=count_succ(v)

frst=(1:(length(v)-1))
nxt=(2:length(v))
mean((v[nxt])[w[frst]])

L’investigateur n’a même pas besoin de compter précisément les patients. Il lui suffit de se souvenir des derniers patients pour anticiper le prochain. Ce problème existe même en cas de bloc de taille aléatoire, même si le calcul de probabilité n’est plus exact.
Dans les essais multicentriques, avec randomisation centralisée, le problème est réduit par le fait que plusieurs investigateurs peuvent inclure sans communiquer. Même si un investigateur peut avoir inclus 4 patients d’affilée dans le même bras, il est possible que d’autres investigateurs, entre temps, aient inclus des patients, de telle sorte que son information sur la séquence aléatoire est brouillée. Néanmoins, il est fréquent de stratifier la randomisation sur le centre, créant ainsi une liste différente pour chaque centre et laissant donc à l’investigateur autant d’information que dans un essai monocentrique.
La seule randomisation qui garantisse des allocations indépendantes, c’est la randomisation simple !
Autrement, l’information disponible sur le groupe de randomisation dans lequel les premiers patients ont été alloués (accessibles dans un essai en ouvert) permet d’obtenir de l’information sur les patients suivants.

C’est pourquoi, au moins dans les essais ouverts, je conseille d’utiliser la randomisation simple. Le risque de biais est ainsi diminué alors que la puissance n’est presque pas abaissée.

Limites du raisonnement

De même que l’évaluation en aveugle, l’aveugle patient et l’aveugle investigateur, la randomisation n’est qu’un outil de rigueur méthodologique qui n’est pas toujours indispensable. L’article Impact of blinding on estimated treatment effects in randomised clinical trials: meta-epidemiological study ne montre pas une différence majeure (même s’il existe une incertitude non négligeable) d’effet selon la rigueur méthodologie des essais cliniques randomisés. Néanmoins, le cas moyen ne décrit pas le cas particulier. Il me paraît indispensable d’être rigoureux dans les contextes où les investigateurs perdent leur neutralité, comme par exemple, en ce moment (avril 2020) pour l’évaluation thérapeutique de l’hydroxychloroquine dans le traitement du COVID-19. Dans ce cas, la randomisation ouverte est à éviter. La randomisation en double aveugle préservant le secret de la séquence d’allocation, la randomisation par blocs n’est alors pas aberrante.