Pourquoi un modèle multivarié ?

Les modèles de régression logistique multivariés se retrouvent dans beaucoup de publications, sans que la raison de leur réalisation soit clairement annoncée. Cela est dommage parce que la manière de les construire dépend beaucoup de cette raison.

Nous listons ici, un nombre de raisons de réaliser un modèle multivarié :

  1. Analyses de causalité en épidémiologie analytique, que ce soit la recherche d’un effet nocif d’une exposition, ou bénéfique d’un traitement
  2. Optimisation de la puissance dans un essai clinique randomisé
  3. Construction d’un score diagnostique ou de dépistage
  4. Recherche d’un modèle pronostique ou prédictif
  5. Recherche de « facteurs de risque »
  6. Redressement d’un échantillon soumis à des biais de sélection
  7. Standardisation directe pour « comparer » des populations entre elles
  8. Imputation simple ou multiple
  9. Interprétation conditionnelle à des variables qu’on ne peut ignorer

1. Analyse de causalité

Dans ces analyses, on cherche à identifier un lien causal entre une exposition et un outcome, et à rechercher éventuellement l’intensité de ce lien. Le concept, c’est qu’en imposant ou supprimant volontairement l’exposition, on espère influencer l’outcome. La recherche de causalité n’a généralement pas d’intérêt sur des facteurs non modifiables. Si on prend la définition contra-factuelle de la causalité, le concept même n’est pas applicable à des facteurs non modifiables. Pour d’autres définitions, il reste applicable.

Il faut ajuster sur les facteurs de confusion, qui sont chronologiquement antérieurs à l’exposition et à l’outcome et sont à la fois des causes de l’exposition et de l’outcome. D’une manière générale, il ne faut pas ajuster sur les variables chronologiquement postérieures à l’exposition et/ou l’outcome. Il ne faut surtout pas ajuster sur les symptômes et/ou conséquences de l’outcome. Il ne faut pas non plus ajuster sur les facteurs de médiation, sauf lorsque cette analyse est complètement assumée (analyse de médiation).

Dans ces modèles, on ne s’intéressera qu’à l’effet de l’exposition sur l’outcome, ajusté sur les covariables. On voudra aussi avoir une estimation non biaisée de cet effet, un intervalle de confiance non biaisé et un petit p interprétable. On voudra avoir une estimation interprétable et quantifiée, telle qu’une différence absolue de risque plutôt qu’un risque relatif dont l’importance dépend beaucoup du risque de base. La modélisation des effets des covariables, par contre, pourra être fine et complexe puisque leurs effets ne seront pas directement interprétés. On fournira des effets marginaux plutôt que conditionnels.

Les approches guidées par les données (data driven), sont, d’une manière générale, biaisées dans ce contexte. Par exemple, l’inclusion de variables corrélées à l’outcome avec p<0.20 en analyse « univariée », va engendrer des fluctuations d’échantillonnages chaotiques rendant invalide les petits p et intervalles de confiance, en plus d’ajuster insuffisamment l’effet d’intérêt. Les techniques automatiques stepwise/backward/forward sont tout aussi biaisées, de même toute technique « manuelle » de recherche du « meilleur modèle ». Par ailleurs, ces méthodes sont susceptibles d’induire en erreur le statisticien en lui faisant ajuster l’effet d’intérêt sur des variables de médiation.

2. Optimisation de la puissance dans un essai clinique randomisé

Pour commencer, il n’existe pas de facteur de confusion dans un essai clinique randomisé bien mené. Le déséquilibre aléatoire des facteurs pronostiques entre les groupes est à l’origine d’une erreur aléatoire d’espérance nulle et pas d’un biais. Évitez de parler de « biais de confusion résiduel », c’est un vocabulaire inadapté. En même temps que de choisir les variables d’ajustement lors de la rédaction d’un protocole d’essai clinique randomisé, il faut commencer par se demander si l’ajustement est pertinent, tout court. Pour cela, se poser les questions :

  1. Qu’est-ce que je gagne ?
  2. Qu’est-ce que je perds ?

Qu’est-ce que je gagne : S’il existe des variables pronostiques majeures, on peut réduire jusqu’à environ 10% le nombre de sujets nécessaires. S’il n’existe que des variables pronostiques « mineures » (p.e. stade TNM IIIa vs IIIb vs IIIc dans un essai clinique incluant des patients avec cancer du colon au stade IIIa), on ne réduira le NSN de 5%, voire moins. Le risque principal, c’est le P-hacking si on a pas été extrêmement explicite sur les procédures d’ajustement dans le protocole et qu’on ne s’est pas assuré que le protocole complet a été publié avant de démarrer l’étude. Il existe d’autres risques, comme le fait de s’orienter vers le modèle logistique sans retransformer l’effet final en différence absolue marginale de risque ou en interprétant un effet conditionnel comme s’il était marginal, ou s’orienter vers un modèle identité-binomial qui est asymptotiquement correct mais sur des échantillons de taille petite ou moyenne est fortement biaisé, d’autant plus qu’on ajuste sur des facteurs fortement pronostiques (cf article de ce blog intitulé « modèle identité-binomial vs identité-gaussien »).

Si on sait ce qu’on fait, alors l’ajustement doit se faire sur des variables suffisamment fortement pronostiques pour que l’erreur d’estimation de l’effet de la variable soit inférieur à l’effet réel de la variable. En bref, le rapport signal/bruit est favorable dans l’estimation de cette variable. Si vous avez une grosse étude, vous pourrez mettre des facteurs faiblement pronostiques. Sur une petite étude, il faut juste mettre les facteurs fortement pronostiques. C’est un pari. Si on ajuste sur des variables peu pertinentes, on augmente un peu l’erreur statistique et baisse un peu la puissance. Si on oublie d’ajuster sur une variable fortement pronostique, on perd aussi un peu en puissance. Les gains et pertes sont minimes (± 1 ou 2% de puissance), sauf pour les variables pronostiques majeures.

Bien sûr, toutes les variables d’ajustement doivent être choisies a priori et les détails méthodologiques hyper-précis doivent être mis dans le protocole, tel que l’estimateur exact utilisé, la méthode d’imputation de chacune des covariables. À la limite, il faut écrire le script d’analyse à l’avance.

Toutes les variables d’ajustement doivent être collectées AVANT la randomisation. Sinon, on risque d’ajuster sur des facteurs de médiation ou de créer des biais d’immortalité et autres variantes.

3. Construction d’un score diagnostique ou de dépistage

La question de la chronologie persiste. Je suppose qu’il existe un Gold Standard, coûteux ou invasif qu’on veut remplacer par des outils diagnostiques plus simples, moins performant ou aussi performant. Dans ce contexte, la question de causalité ne se pose plus. Il n’existe plus de facteur de confusion ou de médiation. Les notions apparaissant sont celles de l’information, de redondance, de « coût » (en temps, en unités monétaires, en effets secondaires, en acceptabilité pour le soignant et patient) de collecte de la variable, mais aussi de parcimonie et de sur-entraînement.

L’usage de méthodes automatiques d’estimation des coefficients, telles que la régression LASSO, ElasticNet ou le stepwise/backward/forward elimination, est autorisée mais n’est pas toujours capable de prendre en compte le coût d’une variable. Il faudra un échantillon de validation distinct de l’échantillon d’entraînement, ou une cross-validation ou autre technique de rectification du sur-entraînement. On a le droit d’utiliser des modèles de machine learning très complexes (p.e. réseaux de neurones), même si je suis convaincu qu’on gagne rarement beaucoup en performance diagnostique avec ces techniques lorsque l’information en entrée est pauvre et qu’on perd beaucoup en transparence et en simplicité d’usage.

La notion de petit p et d’intervalle de confiance disparaît. On doit évaluer les performances diagnostiques du modèle dans son ensemble, par des statistiques telles que l’aire sous la courbe ROC. On doit aussi évaluer la calibration du modèle. On fournira une formule incluant l’ordonnée à l’origine (intercept) du modèle.

La chronologie garde toujours de l’importance. On ne peut pas inclure dans un modèle diagnostique, une variable collectée après le passage du Gold Standard. Par exemple, on ne va pas utiliser la réponse à la radiothérapie contre un cancer du colon comme élément du diagnostic du cancer du colon.

4. Recherche d’un modèle pronostique ou prédictif

Le pronostic, c’est la prédiction de l’évolution favorable (amélioration, guérison) ou défavorable (aggravation, complication, décès) d’une maladie déjà établie.

Le scenario est statistiquement assez similaire à celui d’une étude diagnostique, à ceci près que les variables précèdent chronologiquement l’outcome. Les grands principes restent les mêmes. Il ne faut pas non plus parler de facteur de confusion puisqu’il n’est pas question de causalité, sauf dans de rares cas où on s’intéresse à un facteur pronostique modifiable que l’on pense être causal… Auquel cas, on quitte le champ du pronostic pour rentrer dans le champ de la causalité et les principes énoncés plus haut s’appliquent.

Vous aurez compris que quand on s’intéresse aux biomarqueurs, souvent utilisés dans l’évaluation pronostique ou dans la prédiction de réponse au traitement, on supposent nullement qu’ils sont la cause de l’évolution. Autrement, on rechercherait des traitements qui les ciblent directement. Ils ne sont généralement qu’un symptôme d’un état interne.

Je déconseille l’usage des modèles de Cox avec variable dépendante du temps car ceux-ci ne permettent que la prédiction du présent : dernière mesure de la variable dépendante corrélée avec l’outcome.

Bien sûr, vous l’aurez compris, la chronologie est fondamentale. Rien ne sert de prédire le présent ou le passé. Il faut se placer à un moment chronologiquement bien défini (p.e. diagnostic d’une maladie) et utiliser des variables que l’on peut collecter à ce moment là pour prédire l’événement à un horizon temporel bien défini.

Le moment de l’évaluation doit être explicitement pensé. Par exemple, dans le cas d’un infarctus du myocarde, on peut faire une évaluation pronostique au moment du diagnostic de l’IDM par un gold standard (p.e. coronarographie), on peut faire une seconde évaluation au moment de la sortie d’hospitalisation. Les sujets décédés à l’hôpital ne sont pas concernés du tout par la seconde évaluation pronostique alors qu’au contraire ils peuvent être l’objet de la première ! L’information disponible à la sortie d’hospitalisation est riche et moins volatile (p.e. résultat d’une échographie trans-thoracique à la sortie, degré d’insuffisance cardiaque) que l’information au diagnostic, permettant une évaluation pronostique à moyen et long terme bien plus fine qu’il n’aurait été imaginable de faire au moment du diagnostic.

Attention au modèle de Cox, il conduit au « relativisme ». Tendance à oublier que le risque « de base » est très important à connaître. On peut rechercher les facteurs pronostiques de la sclérose latérale amyotrophique (https://doi.org/10.1007/BF00839964) dont le pronostic global est effroyable (médiane de survie de 2,5 ans) mais considérer que les patients de moins de 65 ans ont un « bon pronostic » parce que leur médiane de survie est de presque 3,5 ans, en oubliant que ce dernier chiffre reste bien petit. Le modèle de Cox fait totalement disparaître de l’équation la courbe de survie de base, trompant donc facilement son monde.

En bref, avant d’expliquer la variance d’un facteur, il faut commencer par estimer son espérance.

Comme pour un modèle diagnostique, un modèle pronostique ou prédictif est basé sur la construction d’une formule permettant de calculer un risque. On peut par exemple, fournir une formule calculant l’espérance de vie restante d’un patient ayant une sclérose latérale amyotrophique.

Ça, c’est la théorie, en pratique, j’ai développé mon propre point de vue, issu de mon expérience qualitative des interactions avec les cliniciens. Il vaut ce qu’il vaut, mais que je vais le partager (c’est l’objet de ce blog).

Les formules, pour la majorité d’entre elles ne seront pas utilisées, et les quelques unes qui le seront (p.e. score Apache) serviront plus à la recherche qu’à la pratique clinique de routine, même si ça dépend des cliniciens, certains étant plus friands que d’autres d’outils qui les guident.

Ces formules ont souvent une très faible validité externe, notamment sur le risque de base, parce qu’elles sont faites sur des échantillons non représentatifs de la population cible, avec en plus de nombreux biais de mesures.

Ensuite, les formules issues de modèles multivariés donnent des poids différents à chaque facteur. On s’aperçoit souvent que cette nuance n’a pas grande importance et qu’en comptant simplement le nombre de facteurs présents, on a des performances pronostiques ou prédictives presque identiques. Au mieux, on peut avoir à distinguer les facteurs majeurs (poids 2) des facteurs mineurs (poids 1).

Ensuite, les cliniciens ont une mémoire énorme, mais limitée. Les consultations sont de courte durée et ils doivent rapidement calculer de tête tous les risques. Cela veut dire qu’ils ne vont pas utiliser la calculatrice que vous leur fournirez, mais ce contenteront de compter consciemment ou inconsciemment le nombre de facteurs de bon ou mauvais pronostic du patient. C’est un peu différent pour la prédiction de réponse au traitement dont la gestion est beaucoup plus codifiée, et nécessite d’être validée par des essais cliniques randomisés.

En bref, de votre article, les cliniciens mémoriseront et utiliseront seulement : les facteurs de risque sont A, B et C. On comprend alors, qu’il est fondamental d’inclure toutes les variables qu’ils utilisent déjà dans la pratique dans votre modèle afin de fournir des nouveaux facteurs pronostiques « indépendants » de ceux qui sont déjà connus. Par exemple, l’âge est un facteur de moralité évident. Si vous analysez le pronostic de la démence à corps de Lewy, vous êtes obligés d’ajuster sur l’âge car l’analyse du clinicien sera toujours conditionnelle à l’âge.

Il ne faut pas non plus inclure de variables qui ne seront pas collectées dans la pratique dans vos formules. Une variable un peu difficile à collecter, si elle apporte beaucoup, sera intéressante à mettre dans le modèle, dans l’idée que même si elle n’est pas encore utilisée en routine, elle le deviendra tellement elle est pertinente. Par contre, n’intégrez pas des variables dans votre modèle si vous savez qu’elles ont été collectées pour la recherche mais ne seront pas utilisées dans la routine clinique. En effet, ces variables risquent d’atténuer, voire de faire disparaître les effets d’autres variables, qui elles seront utilisées dans la routine. Autre formulation : les effets des variables dans le modèle, étant tous conditionnels aux autres, ils ne peuvent être interprétés que lorsque toutes les variables sont mesurées.

In fine, les cliniciens ont besoin d’une simple liste de variables qu’ils sauront employer ensuite. Fournissez leur, autant que possible.

5. Recherche de « facteurs de risque »

Beaucoup d’articles parlent de « rechercher des facteurs de risque » sans préciser ce qu’ils entendent par là. Il y a une ambiguïté. Le terme facteur de risque peut avoir une connotation causale ou pas. Si on souhaite potentiellement modifier ces facteurs de risque pour réduire l’apparition de la maladie, alors les principes de l’analyse causale précédemment décrits s’appliquent, sinon on rentre plutôt dans le champ de la prédiction, avec les principes sus-cités qui s’appliquent.

Plaçons nous maintenant dans un contexte de recherche de facteurs de risque sans que la causalité nous intéresse mais seulement la prédiction.

À quoi cela sert-il ?

À un niveau collectif, cela permet d’identifier des groupes à « haut risque » et à « bas risque » afin de n’appliquer des actions de prévention secondaire (dépistage) qu’à certains sous-groupes. C’est pour cela que le dépistage de certains cancers ne se fait pas en-dessous d’un certain âge. Le risque est trop faible en-dessous de cet âge, pour que les contraintes, risques et coûts associés au dépistage soient rentables. Il faut aussi faire attention à ne pas tomber dans le relativisme. On peut souvent s’apercevoir qu’une action de santé publique n’est pertinente pour personne, ou au contraire est pertinente pour tout le monde. C’est le risque absolu du sous-groupe qu’il faut regarder plutôt que le risque relatif du sous-groupe.

On peut aussi renforcer certaines actions dans des sous-groupes. Il me paraîtrait envisageable d’essayer d’augmenter la couverture vaccinale contre l’hépatite B virale des usagers de drogues intraveineuses, en s’aidant éventuellement du réseau des pharmaciens distribuant les seringues à usage unique (ce n’est qu’une idée). Il est à noter que même si c’est un facteur de risque causal, il n’est modifiable que dans une mesure limitée, et donc, il serait interprété dans cette stratégie comme facteur non modifiable, sans que la causalité ne revête d’intérêt.

À un niveau individuel, pour un clinicien, cela sert parfois à orienter le diagnostic. Une toux, une dyspnée et un amaigrissement progressif chez un patient fumeur orienteront plus rapidement le diagnostic vers le cancer du poumon que chez un non fumeur. Il est à noter que le facteur de risque sera plutôt utilisé comme facteur diagnostic et sa mesure sera synchrone, ou presque, du diagnostic. Le facteur de risque, au sens strict, est utilisé ici comme facteur diagnostic.

6. Redressement d’un échantillon soumis à des biais de sélection

On calcule des pourcentages ou moyennes dans un certain nombre de cases (p.e. combinaison d’une tranche d’âge, un sexe, une catégorie socio-professionnelle), puis on en calcule la moyenne pondérée par le taux de sujets dans ces cases selon une population de référence. Cela correspond au processus de standardisation directe, même si je préfère le terme de redressement lorsqu’il s’agit d’effacer un biais de sélection. Si on veut redresser sur beaucoup de variables, on s’aperçoit rapidement qu’on a zéro ou un sujet pour beaucoup de cases. Cela pose des problèmes d’estimation. Une solution consiste à lisser les estimations en faisant des hypothèses d’additivité de linéarité ou de continuité. On utilisera les « prédictions » du modèle dans chaque case plutôt que le pourcentage ou la moyenne empirique de la case. Il paraît plausible de considérer que le risque de décéder du COVID-19 d’un sujet symptomatique de 71 ans et compris entre celui d’un sujet de 70 ans et celui de 72 ans. Avec un effet polynomial ou une spline, on pourra être même assez fin et estimer avec 4 ou 5 paramètres le risque dans une centaine de cases. Si on ajoute au modèle la présence de maladie respiratoire chronique, en faisant l’hypothèse d’absence d’interaction ou en modélisant de manière simplifiée l’interaction, on pourra seulement ajouter une ou deux paramètres pour doubler le nombre de cases modélisées.

Ces modèles sont alors réalisés avec les idées suivantes :

  1. Le modèle sous-jacent peut-être complexe car il n’est pas communiqué in fine. On s’autorise l’usage de splines, polynômes, interactions, fonctions de lien complexes.
  2. Le choix des variables de modélisation s’impose naturellement : il s’agit de l’ensemble des variables connues dans la population de référence sur lesquelles on craint qu’il existe des biais de sélection
  3. Chaque coefficient supplémentaire intégré dans le modèle réduit le biais d’estimation mais augmente l’erreur. Plus l’échantillon est grand plus les erreurs sont petites alors que les biais sont inchangés. On mettra dans le modèle un nombre de coefficients d’autant plus grand qu’on a un grand échantillon.
  4. Toutes les hypothèses de modélisation sont fausses. Rien ne sert de les tester. Les hypothèses sont extrêmement nombreuses et les tests d’écart aux hypothèses peu puissants, conduisant à un problème de multiplicité des tests en situation de sous-puissance, associé à des fluctuations d’échantillonnages chaotiques associées aux changements de modélisation induits. Si vous êtes tentés de tester une interaction, mettez plutôt le terme correspondant. Au pire, ça augmentera plus l’erreur que ça ne diminuera le biais.

7. Standardisation directe pour « comparer » des populations entre elles

La standardisation directe correspond à la même procédure statistique que le redressement d’échantillon pour biais de sélection. Mais là, je fais référence à des standardisations pour lesquels l’échantillon sélectionné est représentatif de la population dont il est issu mais que l’on réalise afin de comparer les résultats (pourcentages, moyennes) à une ou plusieurs autres populations de structure différente (variables démographiques ± variables sociales ± variables médicales).

L’interprétation des résultats est fortement différente de celle d’un redressement pour biais de sélection. Dans le cadre d’un redressement pour biais de sélection, le résultat brut n’est pas extrapolable à qui que ce soit alors que dans le cadre d’une standardisation, le résultat brut est tout à fait extrapolable à la population dont il est issu et au contraire, le résultat standardisé ne l’est plus. La standardisation sert à effacer l’effet des variables de standardisation afin de comparer les différences qui n’y sont pas dues.

Comme exemple, considérons une nation qui constate une prévalence élevée de démence dans sa population. Plus élevée, par exemple que dans la nation voisine. Après avoir éliminé une différence due à un biais de mesure (sous-diagnostic ou définition diagnostique différente), on peut se demander : la différence est-elle due à une population plus âgée ? En standardisant sur l’âge d’une même population de référence (l’une ou l’autre des nations ou la population de référence OMS), si la différence disparaît, alors l’explication est bien celle-là. Sinon, il y a d’autres différences qui justifient l’écart.

Le chiffre brut, non ajusté, reste le vrai chiffre, impliquant une nécessité d’organisation sociale et médicale pour subvenir aux besoins de cette population de patients déments. Pire encore, la variable de standardisation (âge ± sexe) n’est pas (ou peu) maîtrisable, de telle sorte qu’on va devoir vivre avec ce vrai handicap par rapport au pays voisin. Si, au contraire, on constate qu’une différence avec le pays voisin est explicable par une différence dans la répartition d’un facteur de risque modifiable (p.e. alcoolisme chronique), alors on a un espoir (ce n’est pas gagné pour l’alcool) d’amélioration avec une politique de prévention primaire ou secondaire.

L’interprétation se fait donc, en comparant successivement les chiffres successivement :

  1. Chiffres bruts
  2. Chiffres standardisés sur une ou deux variables (p.e. âge et sexe)
  3. Chiffres standardisés sur encore plus de variables (p.e. prévalence de l’alcoolisme chronique)

On décompose ainsi la différence en ses constituants. La vraie différence, reste la différence brute.

8. Imputation simple ou multiple

Là, on utilise un modèle pour imputer une donnée une ou plusieurs fois. D’un point de vue théorique, on a le droit d’utiliser toutes les autres covariables pour imputer une covariable. La nécessité de respect de la chronologie, obligatoire dans la plupart des modèles, disparaît ici. La seule limite qu’on doit se mettre sur le nombre de covariables explicatives pour l’imputation, c’est le problème de rapport signal/bruit. Si on a un échantillon un peu petit, on brouille le modèle en rajoutant des variables peu pertinentes dont l’effet sera trop mal estimé, avec une erreur plus grande que l’effet. Dans le pire des cas, on met beaucoup trop de variables peu pertinentes et l’imputation par modèle est équivalente à une imputation par une valeur aléatoire. Le choix de l’intégration d’une covariable dans le modèle dépend de la taille d’échantillon, de la variance de cette covariable et de son effet escompté sur la variable à imputer. On a le droit de faire des modèles complexes non linéaires et avec interactions puisqu’ils n’auront pas à être interprétés. On retrouve le principe de modélisation de l’équilibre entre le biais et l’erreur que l’on avait décrit pour le redressement d’échantillon et la standardisation directe.

Lors d’une imputation multiple dans un essai clinique, il ne faut pas oublier de mettre le bras de traitement comme covariable expliquant l’outcome, sans quoi, les valeurs imputées auraient tendance à sous-estimer la différence entre les deux bras, sauf si on assume ce biais (pénalisation volontaire des données manquantes pour préférer un conservatisme à un libéralisme).

9. Interprétation conditionnelle à des variables qu’on ne peut ignorer

Il vous est peut-être arrivé de demander à un clinicien quel était, selon lui, le risque approximatif d’un événement (p.e. complication post-opératoire ou survie sans progression) et qu’il réponde « je ne peux pas dire, ça dépend complètement de … ». Pourquoi cette réponse ? Cela arrive lorsqu’un ou plusieurs facteurs sont hautement pronostiques. Selon le sous-groupe, l’évolution est tellement différente que le clinicien ne moyenne jamais leur pronostic. Une caractéristique de la maladie fortement pronostique est elle qu’on pourrait dire que les patients ayant des valeurs différentes de ce facteur pronostique ont des maladies différentes ! Cela suggère alors que l’interprétation et toujours conditionnelle à ce ou ces facteurs. Il faudrait idéalement faire des analyses en sous-groupes, sauf si on n’a pas la puissance et qu’on a de bonnes raisons de penser que l’interaction est faible ou négligeable. Auquel cas, il faut faire un modèle ajusté sur cette covariable et il est pertinent de fournir un effet conditionnel plutôt que marginal.

Au final, cela se trouve souvent dans les critères d’inclusion. On mélange rarement des patients aux caractéristiques trop différentes. Le fait de n’inclure qu’un sous-groupe de patients est une manière d’à la fois ajuster et prendre en compte les interactions sur la variable définissant se sous-groupe.

There will be legacy

As of April 14, 2020, there are 694 programming languages listed on the Wikipedia list of programming languages. This lists programming languages, not implementations and standards, so that C++ is a single programming language, be it GNU C++, Microsoft’s Visual C++, Clang, Intel C++ compiler or C++98, C++03, C++11, C++14, C++17.

Can a legacy programming language die?

Some primitive Microsoft Excel macros (it was not listed on Wikipedia page of programming languages… I just added it) were available from the first version (1.0 in year 1985) to the fourth version (4.0 in year 1992); they were superseded by Visual Basic for Applications (VBA) macros in Excel 5.0 (year 1993), at a time where personal computers were far less available than today. Microsoft failed to kill these old legacy macros 27 years later. In 2020, the latest Office 365 version still supports Excel 4.0 macro and dialog boxes with their security issues and even has kept the old macro and dialog editor! Why? Because, there are probably millions of lines of code running that old legacy language.

Microsoft tried to kill Visual Basic classic when VB.NET was introduced in 2001. At the same time, it stopped to add any feature to VBA. For instance, the mouse wheel support was never added although the mouse wheel appeared in 1998. When converting the code base to 64 bits, Microsoft did not even take the time to convert the set of common controls and common dialogs to 64 bits, meaning that there are fewer features in newest Office 64 bits versions than in older. Also, that is not official, it is obvious that Microsoft tries to kill VBA. But, the actual successor is very recent. Microsoft added JavaScript/HTML Widgets and is going to add Office Scripts which look very much like VBA macros but with a JavaScript programming language. In my opinion, this is a bad move.

Killing a widely used language is almost impossible. COBOL is 61 years old and still there with the latest standard published in 2014. There are billions lines of codes running out and it is one of the most commonly used language with VB, JavaScript and Java.

Now, Microsoft has to maintain three different macro languages: Excel 1-4 macros, VBA and Office Scripts. Excel 1-4 macros may die, or not, in some near future. VBA are there for much more time than Excel 1-4 macros that did not die in 27 years. If Microsoft survives long enough, I expect VBA to live at least 50 years from 2020, but that may be much longer.

A much better move, in my opinion would be to modernize VBA. Add first-class hashes (Collections are too primitive), closures, real exception handling. Get rid of all assembly code that made transition to 64 bits so hard and rewrite that in C++. Add new controls. Because, anyway, you will have to maintain that almost forewer.

Having many programming languages, create as many problems:

  1. Interoperability between programming languages is poor
  2. A programmer must learn many of them and live with all their flaws
  3. The development ressources allocated to developping implementations is proportional to the number of languages to maintain
  4. Same thing for packages/libraries.

Some counter-argument is that each programming language answers to specific needs. Programming languages are tools. The right tool for the right job must be used. That is ridiculous because almost all programming languages are general purpose and they make pretty much the same thing as others. Okay, there are four families that answers to different needs:

Very low level: assembly languages

Low level: C and C++

Intermediate with strong type system: C#, Java

Scripting languages: PHP, Perl, Python, Ruby, JavaScript, etc.

General purpose scripting languages have pretty all the same features. Multiparadigmatic, weakly typed, supporting closures and object oriented programming.

Specialized programming languages are another beast that have its own problem. I will not analyze them there.

So, if you want to create the next new general purpose programming language, think about the legacy you will create and think about putting the same resources at improving an existing language or its libraries.

If the new programming languages introduces a new paradigm that changes completely the way to program, as OOP and functionnal programming did, then, you certainly should create it. If you just expect to fixes minor flaws of other programming languages, you will see that you just add new flaws and that the old flaws will live 50 years as the old programming languages continue their lives.

Test de normalité

L’usage de tests de normalité n’est pas pertinente dans la recherche biomédicale. J’ai un point de vue très tranché sur la question car cet usage me paraît être une aberration à plusieurs sens, que je détaillerai ci-dessous.

1er argument : les lois du vivant ne sont pas normales

Pour rappel, la loi normale est une loi de distribution de variables quantitatives continues, pouvant prendre toute valeur réelle imaginable de moins l’infini (exclu) à plus l’infini (exclu).

Par principe, toute variable toujours positive ne peut pas suivre une loi normale. Une distribution discrète ne peut pas être normale, ni une distribution bornée. Voilà, j’ai déjà fait le tour de presque toutes les variables quantitatives que j’ai pu voir en recherche biomédicale, allant de l’âge aux concentrations moléculaires (toujours positifs), en passant par l’échelle visuelle analogique (EVA) de douleur, les échelles de qualité de vie (discrètes et bornées).

La normalité suppose aussi :

  1. Que la loi est parfaitement symétrique
  2. Que le kurtosis de la loi est exactement égal à 3
  3. Que tous les mouvements d’ordre supérieur sont parfaitement égaux à ceux d’une loi normale
  4. Qu’après centrage et réduction, la loi de distribution soit parfaitement superposable à celui d’une loi normale centrée-réduite

Si vous avez un doute, vous pouvez toujours consulter cet algorithme.

2ème argument : accepter l’hypothèse nulle est invalide

Pour faire bref, un test de normalité est aussi pertinent qu’un test d’immortalité construit ainsi

Hypothèse nulle H0 : le patient est immortel

Hypothèse nulle H1 : le patient est mortel

Procédure : attendre 1 an. Si le patient est décédé dans l’intervalle, rejeter l’hypothèse nulle H0, sinon, accepter le fait que le patient est immortel.

Cela est aberrant sur deux points : d’abord, cela nie le principe du raisonnement par l’absurde fondant les tests d’hypothèse. Il est pertinent de rejeter l’hypothèse nulle en cas de résultat significatif, mais il est erroné de rejeter l’hypothèse alternative (patient mortel) s’il n’est pas décédé sur une année. Un résultat non significatif, ne permet ni de rejeter, ni d’accepter l’une ou l’autre des hypothèses nulle et alternative. Il ne permet rien de conclure. Ensuite, c’est aberrant parce qu’on sait, dès le départ que tous les patients sont mortels (probabilité a priori = 100%). Le raisonnement réalisé lors de l’acceptation de l’hypothèse nulle est pourtant : jusqu’à preuve du contraire (test significatif), tout patient est immortel. Eh bien, j’inverse la charge de la preuve. Si on veut prouver qu’un patient est réellement immortel, il faut arriver à monter un protocole expérimental rejetant l’hypothèse selon laquelle il est mortel.

Il existe encore pire : la réalisation d’un test de normalité sur une variable ne pouvant, par principe, être normale, telle qu’une variable strictement positive (p.e. hémoglobine) ou discrète (p.e. échelle numérique de douleur de 0 à 10). Cela est équivalent à faire l’autopsie d’un patient dans l’espoir d’échouer à prouver qu’il est mortel pour pouvoir ensuite accepter le fait qu’il est immortel.

3ème argument : le test ne reflète pas le biais des analyses paramétriques

Les test de normalité sont souvent utilisés pour guider l’usage de tests paramétriques (p.e. test de Student) ou non-paramétriques (test de Mann-Whitney). Cela est basé sur l’idée que le test de Student repose sur une hypothèse de normalité. Même si ce test a été développé sur cette théorie, il est asymptotiquement non biaisé si les observations sont indépendantes et identiquement distribuées. Ainsi, il suffit que l’approximation normale de la différence de moyennes soit bonne pour que le test de Student soit valide.

Mon propos est le suivant : les tests de normalité mesurent très mal la qualité de l’approximation normale de la différence de moyennes.

En considérant la distribution comme fixe, plus l’échantillon est grand, plus l’approximation normale du Student est bonne, mais plus le test de normalité est puissant et donc tendra à « rejeter » l’usage du test de Student. Donc, pour une distribution donnée, en faisant varier N, on s’aperçoit que le test de normalité oriente vers la décision opposée à la décision pertinente.

Ensuite, pour un N fixé, les biais associés à l’inférence statistique, telle que l’inflation du risque alpha, sont dépendants de la distribution, mais d’une manière qui est mal décrite par le test de normalité.

Considérons, par exemple, le cas d’un test de Student sur séries appariées de supériorité (risque alpha unilatéral à 2.5%) sur un échantillon de taille 6. Si la distribution des différences appariées suit une loi exponentielle (à un décalage près pour que l’espérance soit zéro), alors le test de Shapiro a une puissance de 22% mais le risque alpha unilatéral du Student monte à 11%, soit une multiplication par un facteur supérieur à 4.4. Si la distribution suit une loi Beta(0.1, 0.1) alors la puissance du test est de 77% bien que le risque alpha unilatéral soit seulement multiplié par 1.03 (soit 2.57%, incertitude 2.566% à 2.579%). Si la distribution suit une loi binomiale de paramètres n=2 et p=0.50, alors la puissance du test est de 44% alors que le risque alpha unilatéral est réduit de 28% relatif (2.5% -> 1.79%). Si la distribution suit une loi de Student à 3 degrés de liberté, alors la puissance du test de normalité est de 11% pour un risque alpha unilatéral du Student descendant à 1.88% (incertitude 1.86 à 1.90%) soit une réduction de 25% relatif. Pour information, les distributions sont présentées ci-dessous :

En bref, sur ces quatre exemples, la plus grande puissance du test de Shapiro porte sur la situation où le test de Student est le moins biaisé.

Le test de Shapiro est sensible aux kurtosis très bas (comme la loi beta) et aux distributions discrètes (comme la loi binomiale) alors que le test de Student résiste très bien à ces deux aspects. Le test de Shapiro est modérément sensible à l’asymétrie de la distribution, mais c’est cette dernière qui biaise le plus le test de Student.

Ensuite, les éléments suivants sont à prendre en compte dans l’évaluation des biais d’un test de Student :

  1. En présence de groupes de taille équilibrée (p.e. 15 vs 15), les asymétries des deux moyennes tendent à se compenser et le test de Student sur séries indépendantes est nettement moins biaisé qu’avec des groupes de taille déséquilibrés (p.e. 15 vs 150).
  2. Les conditions de validité d’un test de Student sur séries appariées reposent sur la distribution des différences appariées et pas la distribution de la variable aléatoire prise isolément dans chaque série de mesures ou dans les deux séries poolées. Par exemple, alors que le poids dans une population obèse, a une distribution asymétrique à droite, les différences de poids après moins avant traitement de l’obésité, suivent une loi plutôt asymétrique à gauche, de forme très différente de la première.

Il faut noter que chaque statistique a ses forces et fragilités. La régression passing Bablok (régression non paramétrique), par exemple, va avoir des problèmes de stabilité en cas de distribution discrète, là où la régression linéaire des moindres carrés y résistera très bien mais aura plus de difficultés à rester stable en présence d’outliers. À ce propos, en présence d’une distribution discrète, le test de normalité risquerait d’orienter vers le passing Bablok, car ce dernier est non paramétrique et donc « résistant aux écart à la normalité ».

4ème argument : les statistiques orientées par les statistiques sont ininterprétables

Faire un test de normalité pour orienter des choix statistiques ultérieurs, c’est méconnaître un des principes fondamentaux de la statistique fréquentiste : une même expérience doit conduire à une même statistique. Par ailleurs, il ne faut pas oublier qu’une statistique répond à une question. Une différence de moyennes peut s’estimer par la méthode de Student ou par du bootstrap paramétrique ou non paramétrique. Ce sont des estimateurs valides et asymptotiquement non biaisés (sauf sur la condition d’homoscédasticité de Student). Le test de Mann-Whitney n’est pas un test de comparaison de moyennes et son risque d’erreur de troisième espèce peut atteindre 100% s’il est interprété comme tel. À côté de cela, l’inflation du risque de première espèce, de 2.5% à 11% unilatéral, d’un Student sur une loi exponentielle avec n=6 est futile.

Quelques néologismes

Si vous suivez ce blog, il peut être intéressant d’apprendre quelques néologismes que je suis susceptible d’utiliser de temps à autre.

Petit péter, petit péteur

Un résultat petit pète, lorsqu’il est statistiquement significatif (petit p < seuil de significativité). Ce terme est neutre et ne présume pas de la réalité de l’hypothèse nulle ou de l’hypothèse alternative, contrairement à la significativité qui a tendance, dans la bouche de non statisticiens, à connoter. Le terme peut se décliner dans divers contextes :

« Ça petit-pète de justesse (p=0.049) »

« Ça petit-pète presque (p=0.06) »

« Ça petit-pète plus fort après ajustement (p=0.001) qu’avant ajustement (p=0.03). »

« Après cette réunion, j’ai compris qu’il voulait petit-péter le plus possible sans se préoccuper de la nature ni de la direction du petit pétage »

Le dernier exemple décrit ce que j’appelle un petit péteur. Le petit péteur est susceptible de faire du P-Hacking, mais il peut utiliser d’autres techniques legit, comme l’enfonçage de portes ouvertes.

Le petit péteur a deux principales formes cliniques :

  1. Le petit péteur orienté, ne cherchant qu’à confirmer ses préjugés, rejetant tout résultat allant à l’encontre de ses préjugés et n’acceptant que ceux qui le confirment. Il veut prouver plutôt que chercher.
  2. Le petit péteur publicatoire atteint du germe Sigapsia publicans, qui ne cherche qu’à maximiser sa publicance

Évidemment, cela est très schématique. La sémiologie du petit péteur est riche et les formes cliniques nombreuses.

Chance alpha

La chance alpha, c’est la probabilité de petit péter conditionnellement à l’hypothèse nulle. C’est exactement la même chose que le risque alpha, mais c’est vu du point de vue du petit péteur. Au lieu de minimiser le risque alpha, il voudra maximiser la chance alpha puisque c’est une chance comme une autre.

Publicance

La publicance est la probabilité inconditionnelle de petit péter, qui est donc à peu près égale à la probabilité de publier.

De manière plus formelle :

Publicance = Proba(petit péter) = Proba(H0)×Proba(petit péter | H0) + Proba(H1)×Proba(petit péter | H1) = Proba(H0)×(chance alpha)+(1-Proba(H0))×puissance

Une méthode statistique augmentera la publicance si elle augmente la puissance, la chance alpha ou les deux. Lorsque que la puissance est faible ou que la probabilité de H0 est élevée, le petit péteur se concentrera sur la maximisation de la chance alpha.

Si on tourne les choses sous un autre angle, on peut dire que la chance alpha est la publicance sous l’hypothèse nulle alors que la puissance est la publicance sous l’hypothèse alternative.

Il faut noter que les méthodes statistiques libérales ont tendance à augmenter la puissance en même temps qu’elles augmentent le risque alpha. Dans ces conditions, je dirais qu’elles augmentent la publicance inconditionnellement à H0 ou H1. Le rapport entre la puissance et le risque alpha diminue généralement dans ces conditions mais le petit péteur n’en aura que faire.

Modèle linéaire ajusté dans un essai clinique randomisé

Encore un billet sur les ajustements dans les essais cliniques randomisés ! Pourquoi s’intéresser à ce sujet ? Parce que plus un modèle et complexe, plus il fait d’hypothèses et plus il est susceptible d’engendrer des biais. Pour un gain de puissance, il paraît difficile d’accepter l’ajout de biais dans un essai clinique randomisé dont le design est pensé pour réduire au maximum les biais. Les seuls biais tolérables sont ceux qui restent infimes, aussi bien en cas de respect des hypothèses faites qu’en cas d’écart à ces hypothèses. Je montrerai par la suite que le modèle identité-gaussien est un bon candidat, même s’il présente quelques imperfections.

Cela me paraît évident, mais je vais répéter que les hypothèses d’absence d’interaction et de linéarité sont toujours fausses, à un degré variable. Par exemple, dans un modèle linéaire généralisé, on a le choix entre plusieurs fonctions de lien (logit, probit, cauchit, cloglog, identité, etc.). Si un des modèles est correct (p.e. linéarité respectée), alors, tous les autres sont faux, puisqu’ils sont fondamentalement contradictoires, sauf dans le cas d’un modèle saturé. Plus vraisemblablement, tous les modèles non saturés sont faux car il n’y a aucune raison qu’une fonction mathématique arbitraire colle parfaitement à la réalité. Néanmoins, les écarts à la linéarité peuvent être modestes ou forts.

Le titre initial de ce billet était « ajustement sur un effet non linéaire dans un essai clinique randomisé ». La théorie initiale de l’auteur, c’était que l’estimateur des moindres carrés du modèle gaussien résistait bien à l’ajustement sur des effets non linéaires, alors que les modèles identité-binomiaux donneraient trop de poids aux sous-groupes où le risque est proche de 0% ou 100% parce que leur variance estimée est trop faible, biaisant les estimations, y compris celles de l’effet du traitement. En pratique, l’estimateur identité-binomial s’est mal comporté même pour des ajustements linéaires. Selon la nature de la non linéarité, les choses peuvent être encore aggravées ou pas. Au total, les estimateurs associés au modèle identité-binomial sont juste plus mauvais que ceux du modèle gaussien, du moins dans le cadre d’un ajustement sur une covariable non corrélée à l’exposition, c’est-à-dire dans le cadre d’essais cliniques randomisés. Au vu de ces résultats, le titre de ce billet a été modifié pour refléter le message.

Nous allons donc nous intéresser à évaluer les biais de différents tests statistiques dans des modèles non ajustés et ajustés sur un facteur pronostique dans le cadre d’essais cliniques randomisés. Un facteur pronostique de forte variance mais de corrélation avec l’outcome variable sera analysé. Les situations de linéarité et de non linéarité seront analysés. Certains effets et écarts à la linéarité seront très forts, de telle sorte que les biais que nous allons décrire seront supérieurs à ceux qu’on rencontrera la plupart du temps. C’est en torturant ainsi les tests que nous identifieront lesquels sont les plus robustes.

Facteur bayésien et puissance corrigée

Généralement, un test qui conduit à un plus gros risque améliore aussi la puissance. Cela est simplement dû au fait que le test est « trop facilement » significatif, peu importe l’hypothèse nulle ou alternative. Afin de pouvoir comparer la puissance deux procédures dont le risque alpha diffère, je définis la notion de « test corrigé », comme un test dont le seuil de significativité nominal est corrigé afin que le risque alpha empirique soit exactement égal à 5%. Le seuil de significativité corrigé est égal au 5ème percentile de la distribution empirique des petits p sous l’hypothèse nulle. Ce test corrigé n’est plus biaisé (risque alpha = 5%). La puissance corrigée est égale à la puissance du test corrigé.

En conditions réelles, il est difficile de créer le test corrigé, car la distribution exacte du facteur pronostique et son lien avec l’outcome ne sont qu’approximativement connus. En pratique, on utilisera le test non corrigé. En conséquence, on vivra avec une inflation du risque alpha et une légère augmentation de la puissance. Par exemple, on passera d’un risque alpha de 5% à 5.1% et d’une puissance de 80% à 85%. Peut-on alors comparer encore les tests qui n’ont pas le même risque alpha ? La réponse est oui, à peu près. Pour ce faire, définissons le facteur bayésien comme le rapport entre la puissance et le risque alpha : FB=Proba(résultat significatif | H1)/Proba(résultat significatif | H0). Cela peut s’interpréter comme un odds ratio de preuve en faveur de H1 lorsque le résultat est significatif. En effet on peut noter Cote(H1)=Proba(H1)/(1-Proba(H1))=Proba(H1)/Proba(H0) la cote de crédibilité a priori de H1, ce qui permet de reformuler le théorème de Bayes comme : Cote(H1 | résultat significatif)=Cote(H1)×FB. Si plusieurs études indépendantes accumulent des preuves, il suffit de multiplier successivement la cote de H1 par les différents facteurs bayésiens afin de calculer la probabilité a posteriori de H1. La proportion de « fausses découvertes » (erreurs de type I) parmi les résultats statistiquement significatifs est fonction de la distribution des probabilités a priori des hypothèses testées (si on cherche n’importe quoi, on trouve n’importe quoi) et des facteurs bayésiens. En considérant que le jeu d’hypothèses testées est figé, un test de meilleur facteur bayésien va fournir une moindre proportion de fausses découvertes.

Cela montre que le facteur bayésien est une statistique pertinente pour la comparaison de tests, mais il a une limite. À facteur bayésien égal, on préfèrera un test plus puissant car, en plus de fournir le même ratio signal/bruit (vraie / fausse découverte), il fournira un plus grand nombre de découvertes, et donc, un plus grand nombre de vraies découvertes en tout. Lorsqu’un test statistique diminue le facteur bayésien mais augmente la puissance, il peut parfois augmenter le nombre de vraies découvertes mais diminuer le ratio vrai/fausse découverte. Dans ces conditions, on ne peut pas dire qu’un test soit complètement supérieur à l’autre.

Exemple : un test avec une puissance à 99% et un risque alpha à 1% est préférable à un test avec une puissance à 99.8% et un risque alpha à 5% car le facteur bayésien est presque cinq fois plus grand, et le nombre de vraie découvertes est à peine diminué.

Par analogie aux tests diagnostiques, le facteur bayésien reflète la VPP/(1-VPP) du test statistique alors que la puissance reflète la sensibilité du test.

Scenarii de simulation

Essai clinique en randomisation simple 1:1 avec N=100, 200 ou 400 patients randomisés

Risque alpha bilatéral nominal : 5%

Outcome binaire : succès vs échec

Ajustement sur un facteur pronostique quantitatif discret, valant 0, 1 ou 2, qui en l’absence de traitement conduit à respectivement, x%, y% et z% de succès sur l’outcome. Plusieurs combinaisons de x, y et z, linéaires et non-linéaires ont été testées.

La distribution du facteur pronostique est : (40%, 40%, 20%) pour les valeurs (0, 1, 2) respectivement.

L’effet du traitement est identique dans chacun des sous-groupes de facteur pronostique.

Deux effets du traitement seront analysés :

  1. Absence totale d’effet (hypothèse « nulle »)
  2. Effet conduisant approximativement à une puissance de 80% pour le modèle gaussien non ajusté (+27.85%, +19.69% ou +13.92% selon que N=100, 200 ou 400) dans chacun des trois sous-groupes pronostiques (hypothèse « alternative »)

Estimateurs comparés :

  1. Modèle identité-gaussien (estimateur : moindres carrés) ajusté sur le facteur pronostique (modèle linéaire général) avec usage d’une statistique t de Student à N-3 degré de liberté pour le test de l’effet traitement
  2. Modèle identité-gaussien non ajusté (équivalent à un test de Student)
  3. Modèle identité-binomial (estimateur : max de vraisemblance) avec test de Wald (sur statistique z approximée à une loi normale centrée réduite) non ajusté
  4. Modèle identité-binomial (estimateur : max de vraisemblance) avec test de Wald (sur statistique z approximée à une loi normale centrée réduite) ajusté sur le facteur pronostique
  5. Modèle identité-binomial (estimateur : max de vraisemblance) avec test du rapport de vraisemblance avec approximation à une loi du chi² de la différence des déviances sous l’hypothèse nulle, non ajusté
  6. Modèle identité-binomial (estimateur : max de vraisemblance) avec test du rapport de vraisemblance avec approximation à une loi du chi² de la différence des déviances sous l’hypothèse nulle, ajusté

Les résultats de tous ces estimateurs ont la même interprétation : différence absolue du risque attribuable au traitement. Par ailleurs, les effets sont supposés se combiner de la même manière (additivité sans transformation) dans ces modèles.

Un modèle identité-binomial et un modèle identité-gaussien avec les mêmes coefficients, prédisent les mêmes espérances aux mêmes observations. Ils ne sont pas tout à fait équivalents car la distribution des observations diffère selon eux. Le modèle identité-gaussien prédit une distribution gaussienne pour la variable à expliquer alors que le modèle identité-binomial prédit une distribution binomiale.

Nous avions déjà montré dans un billet précédent que le modèle identité-gaussien était robuste aux interactions entre le traitement et le facteur pronostique. Nous présumons aussi que le comportement de ce modèle sera bon dans le cas d’un effet non linéaire d’un facteur pronostique car l’estimateur des moindres carrés.

Nombre de simulations : deux millions. Les mêmes données aléatoires ont servi aux 6 tests, de telle sorte que des comparaisons sur séries appariées ont pu être faites pour comparer les performances des estimateurs.

Les incertitudes de simulation, dues au nombre fini de simulations, sont calculées par des intervalles de confiance à 95% selon la méthode de Student à un seul échantillon, à partir de l’échantillon des simulations.

Fréquence de l’outcome

Afin de faciliter la comparaison des résultats, un maximum de scenarii étaient réalisés avec une probabilité marginale de 40% de l’outcome dans le groupe contrôle. Ainsi, on peut considérer que ces scenarii apparaissent dans la même étude, selon le choix du facteur pronostique d’ajustement : faible ou fort, linéaire ou pas.

Gestion des défauts de convergence

Comme la loi binomiale est bornée, certains jeux de coefficients sont impossibles dans un modèle binomial. Par exemple, on ne peut avoir une ordonnée à l’origine à 50%, un facteur pronostique à 1 qui ajoute 30% et un traitement qui ajoute 30%, puisque la combinaison de ces trois effets conduirait à un taux de succès égal à 110%.

Même si le modèle réel est valide, les fluctuations d’échantillonnage conduisent à des fluctuations des coefficients estimés et il arrive qu’il n’y ait plus de jeu de coefficients valide maximisant la vraisemblance de l’échantillon. L’estimateur du maximum de vraisemblance ne converge pas, au sens où la vraisemblance est encore pentue (dérivée non nulle) alors que des observations ont des risques prédits égaux à 100%. Les logiciels statistiques génèrent généralement un message d’erreur dans ces conditions associé ou non au dernier jeu de coefficients avant constatation de la non-convergence.

On peut noter que dans un modèle non ajusté, il n’y a théoriquement toujours un maximum de vraisemblance. Au pire, un des deux groupes atteint 0% ou 100%, mais le maximum de vraisemblance existe toujours. En pratique, les logiciels statistiques ont des difficultés algorithmiques lorsque le maximum de vraisemblance touche ainsi les bords du domaine de définition.

Que ferait le statisticien en charge de l’analyse statistique de l’essai clinique en cas de non-convergence ?

1ère option) Utiliser un estimateur alternatif « de secours », planifié dans le protocole. Cela comprend, par exemple, l’usage du modèle identité-binomial non ajusté, ou l’usage d’un modèle identité-gaussien ajusté ou non, ou encore un test du chi².

2ème option) Utiliser un estimateur alternatif, non planifié dans le protocole car ce cas de figure n’avait pas été anticipé. Il existe alors un risque de P-hacking.

3ème option) Renoncer à fournir un résultat.

Dans le cadre d’essais cliniques randomisés, nous considérons que la 1ère option est la meilleure. Le meilleur modèle de substitution, a priori, est le modèle identité-gaussien avec les mêmes ajustements. C’est ainsi que nous l’avons utilisé dans nos simulations.

Au sens strict, nous n’analysons pas le test du modèle identité-binomial ajusté, mais plutôt la séquence :

  1. Test de Wald ou RV dans un modèle identité-binomial
  2. En cas de défaut de convergence, test de Wald dans un modèle identité-gaussien sur les mêmes données

Valeurs de départ et algorithme

Dans certaines situations, il existe une solution du maximum de vraisemblance tout à fait identifiable, mais le logiciel statistique part de valeurs trop éloignées des coefficients du modèle et se retrouve bloqué au bord du domaine de définition. Pour résoudre ce problème, le choix des coefficients de départ a une grande importance.

Afin d’accélérer la convergence, une première tentative était réalisée avec les coefficients obtenus dans le modèle gaussien avec les mêmes ajustements, via la méthode de décomposition QR à colonne pivot (méthode 0) fournie par la fonction fastglm du package du même nom pour le logiciel R (version 3.6.3). Comme les coefficients du maximum de vraisemblance dans le modèle identité-binomial sont souvent proches de ceux du modèle gaussien, cela permet une convergence rapide. Cependant, cela est aussi à l’origine de défauts de convergence lorsque le modèle gaussien a des prédictions en dehors de l’intervalle [0,1]. En cas d’échec, une deuxième tentative était réalisée avec un intercept à 0.50 et tous les autres coefficients à zéro, prédisant ainsi une chance de succès de 50% pour tous les sujets. Ces coefficients garantissent un éloignement aussi grand que possible des bords du domaine [0,1] et sont associés à une très bonne convergence lorsqu’elle est possible.

Enfin, l’algorithme de fastglm n’est pas toujours optimal. En cas d’échec de cet algorithme avec les deux jeux de valeurs de départ, il était fait appel à la fonction standard stats::glm basée sur la méthode IRWLS (Iteratively ReWeighted Least Squares) avec le deuxième jeu de valeurs de départ prédisant une chance de succès de 50% pour tous les sujets.

Un défaut de convergence était identifié par:

  1. Une erreur lors du fit
  2. Le champ $converged dans l’objet renvoyé par la procédure, valant FALSE
  3. La présence de NA dans les résultats : coefficients du modèle ou petit p.

Résultats modèle Gaussien

N = 100, effet traitement = +27.85 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0510 [0.0507-0.0513]0.0506 [0.0503-0.0509]0.807 [0.806-0.808]0.800 [0.800-0.801]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0510 [0.0507-0.0513]0.0503 [0.0500-0.0506]0.807 [0.806-0.808]0.840 [0.840-0.841]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0510 [0.0507-0.0513]0.0502 [0.0499-0.0506]0.807 [0.806-0.808]0.841 [0.840-0.841]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0510 [0.0507-0.0513]0.0503 [0.0500-0.0506]0.807 [0.806-0.808]0.822 [0.822-0.823]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0510 [0.0507-0.0513]0.0502 [0.0499-0.0505]0.807 [0.806-0.808]0.883 [0.883-0.884]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0510 [0.0507-0.0513]0.0503 [0.0500-0.0506]0.807 [0.806-0.808]0.804 [0.804-0.805]
faible incidence
F=(0.20, 0.20, 0.20)
0.0506 [0.0503-0.0509]0.0499 [0.0496-0.0502]0.847 [0.847-0.848]0.845 [0.844-0.845]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0529 [0.0526-0.0532]0.0511 [0.0508-0.0514]0.830 [0.830-0.831]0.831 [0.830-0.831]
N = 200, effet traitement = +19.69 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0504 [0.0501-0.0508]0.0503 [0.0500-0.0506]0.801 [0.801-0.802]0.797 [0.796-0.797]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0504 [0.0501-0.0508]0.0502 [0.0499-0.0505]0.801 [0.801-0.802]0.833 [0.833-0.834]
effet linéaire fort
F=(0.16, 0.46, 0.76)
0.0504 [0.0501-0.0508]0.0502 [0.0499-0.0505]0.801 [0.801-0.802]0.882 [0.882-0.883]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0504 [0.0501-0.0508]0.0501 [0.0498-0.0504]0.801 [0.801-0.802]0.834 [0.834-0.835]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0504 [0.0501-0.0508]0.0499 [0.0496-0.0502]0.801 [0.801-0.802]0.818 [0.818-0.819]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0504 [0.0501-0.0508]0.0500 [0.0497-0.0503]0.801 [0.801-0.802]0.875 [0.874-0.875]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0504 [0.0501-0.0508]0.0500 [0.0497-0.0503]0.801 [0.801-0.802]0.800 [0.800-0.801]
fort effet en U inversé
F=(0.20, 0.70, 0.20)
0.0504 [0.0501-0.0508]0.0501 [0.0498-0.0504]0.801 [0.801-0.802]0.802 [0.802-0.803]
faible incidence
F=(0.20, 0.20, 0.20)
0.0502 [0.0499-0.0505]0.0499 [0.0496-0.0502]0.867 [0.866-0.867]0.866 [0.865-0.866]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0516 [0.0513-0.0519]0.0505 [0.0502-0.0508]0.812 [0.811-0.812]0.812 [0.811-0.812]
N = 400, effet traitement = +13.92 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0504 [0.0501-0.0507]0.0503 [0.0500-0.0506]0.798 [0.798-0.799]0.797 [0.796-0.797]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0504 [0.0501-0.0507]0.0501 [0.0498-0.0504]0.798 [0.798-0.799]0.834 [0.833-0.834]
effet linéaire fort
F=(0.16, 0.46, 0.76)
0.0504 [0.0501-0.0507]0.0500 [0.0497-0.0503]0.798 [0.798-0.799]0.880 [0.880-0.881]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0504 [0.0501-0.0507]0.0500 [0.0497-0.0503]0.798 [0.798-0.799]0.834 [0.834-0.835]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0504 [0.0501-0.0507]0.0500 [0.0497-0.0503]0.798 [0.798-0.799]0.818 [0.817-0.818]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0504 [0.0501-0.0507]0.0500 [0.0497-0.0503]0.798 [0.798-0.799]0.873 [0.872-0.873]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0504 [0.0501-0.0507]0.0502 [0.0499-0.0505]0.798 [0.798-0.799]0.800 [0.800-0.801]
fort effet en U inversé
F=(0.20, 0.70, 0.20)
0.0504 [0.0501-0.0507]0.0499 [0.0496-0.0502]0.798 [0.798-0.799]0.802 [0.802-0.803]
faible incidence
F=(0.20, 0.20, 0.20)
0.0502 [0.0499-0.0505]0.0501 [0.0498-0.0504]0.884 [0.884-0.885]0.883 [0.883-0.884]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0505 [0.0502-0.0508]0.0502 [0.0499-0.0505]0.803 [0.802-0.803]0.803 [0.802-0.804]

Ces résultats montrent une légère inflation du risque alpha du modèle gaussien non ajusté avec N=100. Cette inflation du risque alpha n’est pas négligeable (0.529) pour une incidence de l’outcome à 50%, mais reste inférieure à +10% relatif (risque alpha < 0.055). Elle est plus faible pour des incidences écartées de 50%. Le modèle ajusté est légèrement plus conservatif tout en étant un peu plus puissant (jusqu’à +8% pour un ajustement sur un fort effet non linéaire) sauf, bien sûr, pour une variable pronostique d’effet nul.

Pour les échantillons plus grands (N=200 et N=400), les effets des traitements ont été réduits afin de maintenir une puissance à 80% environ. Certains scenarii possibles avec N=200 et N=400 n’ont pas été réalisés avec N=100 car certains sous-groupes auraient une incidence de l’outcome au-delà de 100% selon le modèle. Or, nous voulons évaluer les performances des estimateurs, sous l’hypothèse que le modèle est exact, ou invalide seulement sur les défauts de linéarité.

Avec N=200, les risques alpha sont tous très modestes. La puissance gagnée grâce à l’ajustement semble faiblement dépendante de la taille d’échantillon. Ainsi, il semble être aussi intéressant d’ajuster sur un facteur pronostique que l’on ait 400 sujets que 100. C’est compréhensible par le fait que la réduction de variance est d’un facteur assez constant.

Résultats Wald sur identité-binomial

N = 100, effet traitement = +27.85 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0557 [0.0554-0.0560]0.0589 [0.0586-0.0592]0.806 [0.805-0.807]0.795 [0.795-0.796]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0557 [0.0554-0.0560]0.0643 [0.0640-0.0646]0.806 [0.805-0.807]0.779 [0.778-0.779]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0557 [0.0554-0.0560]0.0693 [0.0690-0.0697]0.806 [0.805-0.807]0.771 [0.771-0.772]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0557 [0.0554-0.0560]0.0590 [0.0587-0.0594]0.806 [0.805-0.807]0.824 [0.824-0.825]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0557 [0.0554-0.0560]0.0937 [0.0933-0.0941]0.806 [0.805-0.807]0.737 [0.736-0.738]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0557 [0.0554-0.0560]0.0596 [0.0593-0.0599]0.806 [0.805-0.807]0.793 [0.793-0.794]
faible incidence
F=(0.20, 0.20, 0.20)
0.0559 [0.0556-0.0563]0.0742 [0.0738-0.0745]0.848 [0.848-0.849]0.803 [0.802-0.803]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0566 [0.0562-0.0569]0.0582 [0.0579-0.0585]0.830 [0.829-0.830]0.827 [0.826-0.827]
N = 200, effet traitement = +19.69 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0528 [0.0525-0.0531]0.0539 [0.0536-0.0542]0.801 [0.801-0.802]0.795 [0.795-0.796]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0528 [0.0525-0.0531]0.0553 [0.0550-0.0556]0.801 [0.801-0.802]0.834 [0.833-0.834]
effet linéaire fort
F=(0.16, 0.46, 0.76)
0.0528 [0.0525-0.0531]0.0596 [0.0593-0.0600]0.801 [0.801-0.802]0.825 [0.824-0.825]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0528 [0.0525-0.0531]0.0575 [0.0572-0.0578]0.801 [0.801-0.802]0.822 [0.821-0.823]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0528 [0.0525-0.0531]0.0532 [0.0529-0.0535]0.801 [0.801-0.802]0.825 [0.824-0.825]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0528 [0.0525-0.0531]0.0744 [0.0740-0.0747]0.801 [0.801-0.802]0.825 [0.824-0.825]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0528 [0.0525-0.0531]0.0542 [0.0539-0.0546]0.801 [0.801-0.802]0.797 [0.796-0.797]
fort effet en U inversé
F=(0.20, 0.70, 0.20)
0.0528 [0.0525-0.0531]0.0550 [0.0547-0.0553]0.801 [0.801-0.802]0.797 [0.796-0.797]
faible incidence
F=(0.20, 0.20, 0.20)
0.0527 [0.0524-0.0530]0.0583 [0.0580-0.0587]0.867 [0.867-0.867]0.857 [0.856-0.857]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0539 [0.0536-0.0543]0.0539 [0.0536-0.0542]0.812 [0.811-0.812]0.812 [0.811-0.812]
N = 400, effet traitement = +13.92 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0516 [0.0513-0.0519]0.0520 [0.0517-0.0523]0.798 [0.798-0.799]0.796 [0.796-0.797]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0516 [0.0513-0.0519]0.0523 [0.0520-0.0526]0.798 [0.798-0.799]0.836 [0.835-0.836]
effet linéaire fort
F=(0.16, 0.46, 0.76)
0.0516 [0.0513-0.0519]0.0538 [0.0534-0.0541]0.798 [0.798-0.799]0.894 [0.894-0.894]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0516 [0.0513-0.0519]0.0542 [0.0538-0.0545]0.798 [0.798-0.799]0.831 [0.831-0.832]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0516 [0.0513-0.0519]0.0511 [0.0508-0.0514]0.798 [0.798-0.799]0.820 [0.820-0.821]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0516 [0.0513-0.0519]0.0627 [0.0624-0.0631]0.798 [0.798-0.799]0.863 [0.862-0.863]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0516 [0.0513-0.0519]0.0523 [0.0520-0.0526]0.798 [0.798-0.799]0.800 [0.799-0.800]
fort effet en U inversé
F=(0.20, 0.70, 0.20)
0.0516 [0.0513-0.0519]0.0526 [0.0523-0.0529]0.798 [0.798-0.799]0.801 [0.800-0.801]
faible incidence
F=(0.20, 0.20, 0.20)
0.0513 [0.0510-0.0516]0.0534 [0.0531-0.0537]0.884 [0.884-0.884]0.881 [0.880-0.881]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0512 [0.0509-0.0515]0.0516 [0.0513-0.0519]0.803 [0.802-0.803]0.803 [0.802-0.803]

Le Wald identité-binomial avec N=100 marche très mal dès qu’il y a un ajustement. Même sans ajustement, il y a une petite inflation du risque alpha (> +10% relatif). En présence d’ajustement, le risque alpha peut presque doubler par rapport au risque nominal (effet non linéaire fort, alpha>0.093). La non-linéarité ne pose pas tant problème que la force de prédiction du facteur pronostique (cf effet linéaire fort avec N=200 vs effet en U inversé avec N=200). Paradoxalement, une relation en U inversé va conduire à un biais moins fort (0.0542) qu’un effet linéaire fort (0.0596). Cela est explicable par le fait que l’effet estimé de la covariable sera proche de zéro de telle sorte que le modèle ajusté n’est pas très différent du modèle non ajusté. Il semble que l’inflation du risque alpha arrive surtout quand le modèle prédit à la fois des probas basses et des probas élevées (ce qui est majeur avec le modèle non linéaire fort). La puissance corrigée est plus faible que pour les modèles gaussiens. Avec N=100, l’ajustement fait plutôt perdre en puissance corrigée, alors qu’avec N=200 ou N=400, il fait plutôt gagner.

Avec N=400, le comportement est bien meilleur qu’avec N=100 ou N=200, mais dans les cas extrêmes le risque alpha peut quand même atteindre 6.27% (soit +25% relatif par rapport au risque nominal) avec N=400.

Résultats avec RV sur identité-binomial

N = 100, effet traitement = +27.85 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0525 [0.0521-0.0528]0.0541 [0.0538-0.0544]0.806 [0.806-0.807]0.798 [0.797-0.799]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0525 [0.0521-0.0528]0.0568 [0.0564-0.0571]0.806 [0.806-0.807]0.848 [0.847-0.848]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0525 [0.0521-0.0528]0.0575 [0.0572-0.0578]0.806 [0.806-0.807]0.831 [0.830-0.831]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0525 [0.0521-0.0528]0.0547 [0.0544-0.0550]0.806 [0.806-0.807]0.834 [0.833-0.834]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0525 [0.0521-0.0528]0.0645 [0.0642-0.0649]0.806 [0.806-0.807]0.869 [0.868-0.869]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0525 [0.0521-0.0528]0.0537 [0.0534-0.0541]0.806 [0.806-0.807]0.801 [0.800-0.802]
faible incidence
F=(0.20, 0.20, 0.20)
0.0531 [0.0528-0.0534]0.0620 [0.0617-0.0624]0.844 [0.844-0.845]0.826 [0.826-0.827]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0537 [0.0533-0.0540]0.0541 [0.0538-0.0544]0.831 [0.830-0.831]0.831 [0.830-0.831]
N = 200, effet traitement = +19.69 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0514 [0.0510-0.0517]0.0519 [0.0516-0.0522]0.801 [0.801-0.802]0.796 [0.795-0.797]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0514 [0.0510-0.0517]0.0528 [0.0524-0.0531]0.801 [0.801-0.802]0.835 [0.834-0.835]
effet linéaire fort
F=(0.16, 0.46, 0.76)
0.0514 [0.0510-0.0517]0.0555 [0.0552-0.0558]0.801 [0.801-0.802]0.908 [0.907-0.908]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0514 [0.0510-0.0517]0.0528 [0.0524-0.0531]0.801 [0.801-0.802]0.830 [0.829-0.830]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0514 [0.0510-0.0517]0.0519 [0.0516-0.0522]0.801 [0.801-0.802]0.821 [0.821-0.822]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0514 [0.0510-0.0517]0.0574 [0.0570-0.0577]0.801 [0.801-0.802]0.863 [0.862-0.863]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0514 [0.0510-0.0517]0.0515 [0.0512-0.0518]0.801 [0.801-0.802]0.799 [0.799-0.800]
fort effet en U inversé
F=(0.20, 0.70, 0.20)
0.0514 [0.0510-0.0517]0.0517 [0.0514-0.0520]0.801 [0.801-0.802]0.801 [0.800-0.801]
faible incidence
F=(0.20, 0.20, 0.20)
0.0514 [0.0511-0.0517]0.0544 [0.0541-0.0547]0.866 [0.866-0.867]0.860 [0.860-0.860]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0525 [0.0522-0.0528]0.0520 [0.0517-0.0523]0.811 [0.811-0.812]0.812 [0.811-0.812]
N = 400, effet traitement = +13.92 %
Risque alpha
non ajusté
Risque alpha
ajusté
Puissance corrigée
non ajusté
Puissance corrigée
ajusté
effet nul
F=(0.40, 0.40, 0.40)
0.0509 [0.0506-0.0512]0.0510 [0.0507-0.0513]0.798 [0.797-0.799]0.797 [0.796-0.797]
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.0509 [0.0506-0.0512]0.0512 [0.0509-0.0515]0.798 [0.797-0.799]0.836 [0.835-0.836]
effet linéaire fort
F=(0.16, 0.46, 0.76)
0.0509 [0.0506-0.0512]0.0523 [0.0519-0.0526]0.798 [0.797-0.799]0.894 [0.894-0.895]
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
0.0509 [0.0506-0.0512]0.0512 [0.0509-0.0515]0.798 [0.797-0.799]0.835 [0.835-0.836]
non linéarité alternative
F=(0.33, 0.33, 0.68)
0.0509 [0.0506-0.0512]0.0509 [0.0506-0.0512]0.798 [0.797-0.799]0.818 [0.818-0.819]
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.0509 [0.0506-0.0512]0.0528 [0.0525-0.0532]0.798 [0.797-0.799]0.878 [0.877-0.878]
effet en U inversé
F=(0.24, 0.64, 0.24)
0.0509 [0.0506-0.0512]0.0508 [0.0505-0.0511]0.798 [0.797-0.799]0.801 [0.800-0.801]
fort effet en U inversé
F=(0.20, 0.70, 0.20)
0.0509 [0.0506-0.0512]0.0506 [0.0503-0.0509]0.798 [0.797-0.799]0.803 [0.802-0.804]
faible incidence
F=(0.20, 0.20, 0.20)
0.0508 [0.0505-0.0511]0.0520 [0.0517-0.0523]0.884 [0.883-0.884]0.881 [0.881-0.882]
incidence idéale
F=(0.50, 0.50, 0.50)
0.0508 [0.0505-0.0511]0.0508 [0.0505-0.0511]0.802 [0.802-0.803]0.803 [0.802-0.804]

Pour N=100, le test du rapport de vraisemblance est nettement moins biaisé que le test de Wald dans le modèle identité-binomial. La puissance corrigée est proche de celle du modèle gaussien, mais il y a des petites différences.

Le RV identité-binomial ajusté gagnait un peu en puissance corrigée pour le scenario (N=100) de non linéarité alternative (+1.13%, incertitude +1.11% à 1.15%) et l’effet linéaire moyen (+0.78%, incertitude +0.76% à +0.81%). Le modèle identité-gaussien gagnait en puissance corrigée dans les autres scenarii, jusqu’à +1.72% (incertitude +1.70% à +1.74%) pour le scenario de faible incidence (N=100). Pour le scenario d’incidence idéale (F=(50,50,50)), les deux millions de simulations n’ont pas suffi à trouver la différence (incertitude −0.14‰ à +1.65‰ pour RV moins gaussien).

Avec N=200, les biais étaient nettement plus raisonnables sauf dans les scenarii extrêmes où l’inflation du risque alpha pouvait légèrement dépasser +10% relatif.

Avec N=400, toutes les inflations de risque alpha étaient en dessous de la limite de +10% relatif.

Facteur bayésien

N = 100, effet traitement = +27.85 %
Facteur
Bayésien
id-gauss non ajusté
Facteur
Bayésien
id-gauss ajusté
Facteur
Bayésien
id-binom-wald ajusté
Facteur
Bayésien
id-binom-RV ajusté
effet nul
F=(0.40, 0.40, 0.40)
15.915.913.914.9
effet linéaire moyen
F=(0.24, 0.44, 0.64)
15.916.712.515.2
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
15.916.711.714.7
non linéarité alternative
F=(0.33, 0.33, 0.68)
15.916.414.315.4
effet non linéaire fort
F=(0.10, 0.60, 0.60)
15.917.68.613.8
effet en U inversé
F=(0.24, 0.64, 0.24)
15.91613.715.1
faible incidence
F=(0.20, 0.20, 0.20)
16.816.911.513.7
incidence idéale
F=(0.50, 0.50, 0.50)
15.816.314.515.5
N = 200, effet traitement = +19.69 %
Facteur
Bayésien
id-gauss non ajusté
Facteur
Bayésien
id-gauss ajusté
Facteur
Bayésien
id-binom-wald ajusté
Facteur
Bayésien
id-binom-RV ajusté
effet nul
F=(0.40, 0.40, 0.40)
15.915.914.915.4
effet linéaire moyen
F=(0.24, 0.44, 0.64)
15.916.615.315.9
effet linéaire fort
F=(0.16, 0.46, 0.76)
15.917.61416.5
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
15.916.714.615.8
non linéarité alternative
F=(0.33, 0.33, 0.68)
15.916.415.615.9
effet non linéaire fort
F=(0.10, 0.60, 0.60)
15.917.511.715.3
effet en U inversé
F=(0.24, 0.64, 0.24)
15.91614.915.6
fort effet en U inversé
F=(0.20, 0.70, 0.20)
15.91614.715.6
faible incidence
F=(0.20, 0.20, 0.20)
17.317.314.916
incidence idéale
F=(0.50, 0.50, 0.50)
15.816.115.215.7
N = 400, effet traitement = +13.92 %
Facteur
Bayésien
id-gauss non ajusté
Facteur
Bayésien
id-gauss ajusté
Facteur
Bayésien
id-binom-wald ajusté
Facteur
Bayésien
id-binom-RV ajusté
effet nul
F=(0.40, 0.40, 0.40)
15.915.915.415.7
effet linéaire moyen
F=(0.24, 0.44, 0.64)
15.916.616.116.4
effet linéaire fort
F=(0.16, 0.46, 0.76)
15.917.616.717.2
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
15.916.715.516.4
non linéarité alternative
F=(0.33, 0.33, 0.68)
15.916.416.116.1
effet non linéaire fort
F=(0.10, 0.60, 0.60)
15.917.414.116.7
effet en U inversé
F=(0.24, 0.64, 0.24)
15.91615.415.8
fort effet en U inversé
F=(0.20, 0.70, 0.20)
15.916.115.315.9
faible incidence
F=(0.20, 0.20, 0.20)
17.617.616.617
incidence idéale
F=(0.50, 0.50, 0.50)
15.91615.715.9

Comme attendu, le modèle identité-gaussien ajusté avait un facteur Bayésien constamment supérieur aux autres, sauf en l’absence totale d’effet, scenario dans lequel le modèle non ajusté était très légèrement supérieur. Dans le pire des scenarii d’effet non linéaire fort avec N=100, le facteur bayésien de l’identité-gaussien ajusté (17.6) était deux fois plus grand que celui du Wald identité-binomial (8.6). Ce scenario n’est néanmoins pas très réaliste car il est rare qu’on trouve un facteur pronostic aussi fort et aussi éloigné de la linéarité.

Avec un nombre de sujets élevé (N=400), les écarts s’atténuaient beaucoup entre les tests.

Défaut de convergence

N = 100, effet traitement = +27.85 %
Taux de non-convergence sous H0
id-binom-Wald ajusté
Taux de non-convergence sous H0
id-binom-RV ajusté
Taux de non-convergence sous H1
id-binom-Wald ajusté
Taux de non-convergence sous H1
id-binom-RV ajusté
effet nul
F=(0.40, 0.40, 0.40)
3.2e-053.2e-050.000170.00017
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0.000250.000250.0280.028
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
5e-040.000520.0920.092
non linéarité alternative
F=(0.33, 0.33, 0.68)
1e-041e-040.00140.0015
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.00940.00990.390.4
effet en U inversé
F=(0.24, 0.64, 0.24)
1.4e-051.4e-050.00810.0081
faible incidence
F=(0.20, 0.20, 0.20)
0.00310.00310.00290.0029
incidence idéale
F=(0.50, 0.50, 0.50)
3e-063e-060.00180.0018
N = 200, effet traitement = +19.69 %
Taux de non-convergence sous H0
id-binom-Wald ajusté
Taux de non-convergence sous H0
id-binom-RV ajusté
Taux de non-convergence sous H1
id-binom-Wald ajusté
Taux de non-convergence sous H1
id-binom-RV ajusté
effet nul
F=(0.40, 0.40, 0.40)
0000
effet linéaire moyen
F=(0.24, 0.44, 0.64)
000.000890.00089
effet linéaire fort
F=(0.16, 0.46, 0.76)
5.4e-055.5e-050.0280.028
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
1e-061e-060.00480.0048
non linéarité alternative
F=(0.33, 0.33, 0.68)
004.5e-064.5e-06
effet non linéaire fort
F=(0.10, 0.60, 0.60)
0.000320.000340.130.13
effet en U inversé
F=(0.24, 0.64, 0.24)
001.4e-051.4e-05
fort effet en U inversé
F=(0.20, 0.70, 0.20)
006.6e-056.6e-05
faible incidence
F=(0.20, 0.20, 0.20)
9.7e-059.7e-050.000120.00012
incidence idéale
F=(0.50, 0.50, 0.50)
001.5e-061.5e-06
N = 400, effet traitement = +13.92 %
Taux de non-convergence sous H0
id-binom-Wald ajusté
Taux de non-convergence sous H0
id-binom-RV ajusté
Taux de non-convergence sous H1
id-binom-Wald ajusté
Taux de non-convergence sous H1
id-binom-RV ajusté
effet nul
F=(0.40, 0.40, 0.40)
0000
effet linéaire moyen
F=(0.24, 0.44, 0.64)
0000
effet linéaire fort
F=(0.16, 0.46, 0.76)
000.000480.00048
effet non linéaire moyen
F=(0.19, 0.54, 0.54)
002.5e-062.5e-06
non linéarité alternative
F=(0.33, 0.33, 0.68)
0000
effet non linéaire fort
F=(0.10, 0.60, 0.60)
1e-061e-060.00490.0049
effet en U inversé
F=(0.24, 0.64, 0.24)
0000
fort effet en U inversé
F=(0.20, 0.70, 0.20)
0000
faible incidence
F=(0.20, 0.20, 0.20)
0000
incidence idéale
F=(0.50, 0.50, 0.50)
0000

Aucun défaut de convergence n’a été observé pour les modèles gaussiens. Ces modèles convergent dès qu’il y a au moins une observation par groupe de traitement. Le seul problème peut concerner le calcul du test de Wald, avec une division par zéro si la variance résiduelle est nulle, ce qui ne peut arriver qu’au cas où, sur l’échantillon, le taux d’événement est de 0% ou 100%. Dans le scenario le plus défavorable, (F=(0.20, 0.20, 0.20), N=100), la probabilité que cela arrive est égale à dbinom(0,100, 0.20) = 2×10^-10. Les défauts de convergence n’étaient pas rares avec le modèle identité-binomial. Le risque est un peu plus élevé pour le test du RV que pour le test de Wald car ce premier nécessite d’estimer deux fois le modèle (avec et sans l’effet traitement), mais la différence est négligeable. Les non convergences étaient globalement peu fréquentes (< 1%) sauf avec N=100 en cas d’ajustement sur un facteur pronostique avec effet linéaire moyen (2.8%) ou non linéaire moyen (9.2% de risque) ou fort (39 ou 40% de risque). Ces non convergences doivent être anticipées et proprement planifiées dans le protocole, ce qui nécessite une grande rigueur. En cas d’effet non linéaire, le choix du modèle de secours a une importance non négligeable, étant donné qu’il est fréquemment utilisé.

Dans les modèles binomiaux non ajustés, seuls quatre défauts de convergence ont été retrouvés sur les soixante millions de simulations, tous dans le scenario F=(0.20,0.20,0.20). Le modèle gaussien non ajusté a servi comme modèle de secours dans ces cas.

Biais de l’estimation ponctuelle

Avec N=100, 40% d’événements et un effet traitement de +0.2785, l’estimateur des moindres carrés (modèle identité-gaussien) non ajusté semblait fournir une estimation ponctuelle de biais infime sous H1 : le biais était estimé à moins de 10^-16. Avec N=30, ce même estimateur avait un biais égal à 2×10^-9. Avec N=20, le biais montait à 5×10^-7 restant donc absolument négligeable. Ces calculs extrêmement précis ne sont pas issus de simulations mais de calculs « exacts » (aux arrondis FP64 près) basés sur la connaissance de la distribution binomiale (fonction dbinom() du logiciel R). La précision du résultat est d’environ 10^-16.

L’estimateur des moindres carrés ajusté sous H1 fournissait une estimation ponctuelle dont le biais n’a pas pu être mesuré avec cent millions de simulations (< 1.2×10^-5), même en situation d’ajustement sur une variable non linéaire. Afin d’affiner la précision, pour le modèle avec effet non linéaire fort, vingt millions de simulations ont été réalisées. Le modèle non ajusté ne montrait pas de biais d’estimation ponctuelle sous H1 (incertitude -9.6×10^-6 à +7.5×10^-5) et le modèle ajusté n’en montrait pas plus (incertitude -8.3×10^-6 à +3.0×10^-5).

L’estimateur du maximum de vraisemblance sur le modèle identité-binomial semblait peu sujet aux biais en l’absence d’ajustement. En présence d’ajustement avec N=100, il existait une très légère surestimation en présence d’un effet nul (+0.21‰ absolu, incertitude +0.08‰ à +0.34‰) ainsi qu’en présence d’un effet linéaire moyen (+3.26‰ , incertitude +3.13 à +3.39‰). En présence d’effets non linéaires, le biais pouvait changer dans un sens ou l’autre. Le plus gros biais était enregistré pour le scenario de « non linéarité alternative » avec une espérance de l’estimateur augmenté de +4.41‰ (incertitude 4.28‰ à 4.54‰). Ces biais tendaient à s’atténuer avec la taille d’échantillon. En comparaison à l’erreur habituellement retrouvée sur les estimations, ces biais sont infimes.

Discussion

Étant donné que la fréquence des outcomes d’une loi binaire suit une loi binomiale si les observations sont indépendantes, il paraîtrait logique d’estimer les effets dans un modèle binomial. Paradoxalement les estimateurs habituels d’un modèle gaussien (moindres carrés + Wald) est plus performant (moins biaisé, meilleur facteur bayésien) que les estimateurs dédiés aux modèles binomiaux (maximum de vraisemblance + Wald ou RV). Les modèles identité-binomiaux souffrent aussi de problème de non-convergence s’ils sont ajustés sur de forts facteurs pronostiques, obligeant à anticiper un modèle de secours. Si ce dernier n’est pas spécifié dans le protocole, la voie du P-hacking est ouverte.

La simplicité du modèle linéaire gaussien, sa robustesse aux écarts aux conditions de validité (interactions, défaut de linéarité), le faible risque d’erreur de manipulation du logiciel (contrairement aux effets marginaux calculés à partir d’un modèle logistique), en fait un excellent choix. L’ajustement peut faire gagner un peu en puissance et en facteur bayésien, mais ses limites ont déjà été discutées.

Le théorème central limite garantit la convergence asymptotique vers une loi normale des coefficients d’un modèle linéaire gaussien. L’estimateur de la variance peut être biaisé, sans être jamais asymptotiquement correct en cas d’effet non linéaire. C’est pourquoi sur des échantillons de grande taille on peut envisager de bootstrapper pour obtenir des estimations correctes.

Au total, il n’y a pas de raison d’utiliser les modèles identité-binomiaux dans les essais cliniques randomisés. Ils n’ont aucun avantage par rapport au modèle identité-gaussien.

Combien de noeuds pour paralléliser mes simulations ?

Si vous faites des simulations statistiques lourdes, vous devez ou devriez bénéficier de la parallélisation des ordinateurs multi-coeurs. Un quadri-coeur peut exécuter environ quatre fois plus vite (sauf saturation de la bande passante RAM) les simulations qu’un mono-coeur. Les outils de parallélisation tels que le package R « parallel » permettent de choisir le nombre de noeuds d’exécution. Quelle valeur choisir ?

La réponse est simple : choisissez autant de noeuds que votre processeur n’a de threads, ce qui est habituellement égal au double du nombre de coeurs. Vous pouvez éventuellement en mettre un ou deux de moins si vous voulez continuer à utiliser votre PC avec confort.

Exemple : processeur 4 coeurs et 8 threads.

Si vous utilisez quatre noeuds, vous gagnerez en performance d’un facteur quatre.

Si vous utilisez cinq noeuds, vous gagnerez en performance d’un facteur 4,2 environ par rapport à un seul noeud. Et vous gagnerez environ 0,2 unité de performance par noeud supplémentaire jusqu’à 8 noeuds. Avec 8 noeuds vous aurez un calcul 4+0.2×4=4.8 fois plus rapide qu’avec un seul noeud. Au delà, vous perdrez un peu en performance.

Limites de la règle n=threads

Parfois on peut gagner à mettre moins d’unité de calcul que ça, pour ne pas saturer la bande passante de la RAM lorsque les algorithmes travaillent sur de grandes quantité de données à la fois. Il est rare que vous ayez intérêt à en mettre moins que de coeurs.

Si on accède intensément à un disque dur, alors les accès alternés peuvent être très lent, à cause du temps d’accès long, et il peut être préférable de ne lancer qu’un seul noeud. Ce phénomène est moindre avec les Solid State Drives (SSD).

Pourquoi ?

Un coeur est un ensemble d’unités de calcul arithmétique et flottantes regroupés pour exécuter du code. Lorsqu’on utilise un coeur simplement, les unités de calcul ne sont pas saturées 100% du temps car certains calculs dépendent des précédents. Par exemple, si on calcule (5+3)×(2+8+4), on doit réaliser trois additions avant de réaliser la multiplication. Pendant que les unités de calcul d’addition sont utilisées, l’unité de multiplication est susceptible d’être au repos. Si on utilise deux threads sur le même coeur, un deuxième fil d’exécution est susceptible de saisir cette opportunité pour faire une multiplication dont elle a besoin à ce moment là. Néanmoins, avec deux threads sur un même coeur, ce dernier est presque en permanence à saturation et on ne gagne qu’environ 20% par rapport à une exécution séquentielle des simulations sur un seul thread.

Les détails sont plus complexes mais nécessitent de maîtriser l’architecture interne d’un PC. Une unité de calcul est plus ou moins spécialisée (p.e. Unité Arithmétique et Logique ou ALU en anglais).

Faut-il des petits p dans le tableau 1 d’un essai clinique randomisé ?

Introduction

Le tableau 1, dans les essais cliniques randomisés, décrit généralement les caractéristiques initiales des deux groupes randomisés. Âge, sexe, état général, comorbidités, sévérité de la maladie d’intérêt. Les grandes revues, telles que je New England Journal of Medicine présentent généralement les caractéristiques sous forme d’un tableau avec une colonne par groupe de randomisation. Il arrive, dans d’autres revues, que des comparaisons statistiques soient faites entre les deux groupes sous forme d’une dernière colonne de petits p. Quelle présentation est la plus pertinente ?

Résumé

Un petit p entre 0.01 et 0.05 n’a aucune valeur. Un petit p à 10^-9 prouve un écart au protocole de randomisation et est donc fortement informatif. Les petits p sont donc utiles s’ils sont bien interprétés. Afin d’éviter une surinterprétation d’un petit p entre 0.01 et 0.05 par un lecteur mal avisé ou un selective reporting d’auteurs qui souhaiteraient dissimuler un écart au protocole, le compromis du NEJM est optimal : fournir les caractéristiques des patients dans chaque groupe, mais ne pas faire de test statistique formel tout en permettant au lecteur avisé de réaliser ces tests a posteriori.s

Argumentaire détaillé

Dans les études observationnelles, on présente ce tableau avec des petits p afin de « montrer » la « comparabilité » des groupes. Je ne m’étendrai pas trop sur les études observationnelles mais traiterai essentiellement du cas des essais cliniques randomisés.

L’argument pour défendre l’absence de présentations des petits p est le fait que, par principe, toute différence observable entre les deux groupes est due au hasard, cela étant garanti par le processus de randomisation. Ainsi, par exemple, une différence d’âge conduisant à un petit p à 0.01 entre les deux groupes, est entièrement due à des fluctuations d’échantillonnage et même si cela désavantage un des deux groupes, cela n’affecte pas l’espérance de la statistique finale. Cela reflète de l’erreur aléatoire d’espérance nulle, pas du biais. Ceux qui veulent réduire cette erreur aléatoire peuvent préalablement spécifier dans le protocole un ajustement sur des facteurs pronostiques majeurs ainsi qu’une stratification de la randomisation sur ces mêmes facteurs (les limites en ont déjà été discutées dans ce blog). En conséquence, le garant de la comparabilité initiale des groupes est le processus de randomisation et aucunement la non-significativité des petits p à baseline. Donc, les petits p ne servent pas à savoir si les groupes sont comparables à baseline.

Il y a encore un autre problème avec la validité statistique des petits p dans le tableau 1. Dès que la randomisation est stratifiée sur un facteur, l’équilibre du facteur entre les deux groupes est garanti comme presque parfait, de telle sorte que les tests usuels (chi², Student), fournissant des petits p uniformes entre 0 et 1 en cas de randomisation simple, auront tendance à fournir des petits p très proches de 1. D’un point de vue strict, les tests employés sont statistiquement erronés.

Ensuite, le problème de multiplicité des tests existe. Souvent il y a 10 ou 20 caractéristiques présentées, et donc 10 ou 20 tests. Il est donc tout à fait attendu d’avoir un ou deux petits p inférieur à 0.05 dans ce tableau, en conditions tout à fait ordinaires, sans la moindre malchance.

Enfin, aussi bien dans les études observationnelles que randomisées, il y a un problème à réaliser un test de différence lorsqu’on souhaite montrer l’équivalence. Ce problème est caricatural dans les études observationnelles de toute petite taille pour lesquelles des différences énormes sont observées dans les groupes et on conclut que les deux groupes ont les mêmes caractéristiques dans les populations dont elles sont tirées puisque les différences peuvent être expliquées par le hasard. Néanmoins, cela ne prouve pas qu’elles SONT expliquées par le hasard. D’un point de vue fréquentiste, seul un test significatif a valeur de preuve. Échouer à montrer une différence ne prouve pas l’équivalence, notamment lorsque la puissance est presque nulle. Si on veut prouver l’équivalence des caractéristiques des groupes, des tests d’équivalence doivent être réalisés plutôt que des tests de différence.

Néanmoins, c’est là que l’interprétation bayésienne doit intervenir. De mon point de vue, le petit p de différence dans un essai clinique randomisé peut être pertinent à condition de placer le seuil de significativité extrêmement bas (10^-6 pour éveiller le soupçon, 10^-9 pour le confirmer). Il ne faut pas négliger le risque d’écart au protocole de randomisation. La qualité de la randomisation fait partie des éléments à prendre en compte dans la lecture critique d’article, et le tableau 1 peut aider à juger de cette qualité. Ainsi, un investigateur d’un centre où il est le seul à inclure peut aisément tricher dans un essai clinique randomisé en ouvert avec une randomisation par blocs de taille 4 stratifiée sur le centre. En effet, l’aléatoire devient très prévisible. Dès que 4 patients ont été randomisés dans le même bras, l’investigateur sait avec certitude que le patient suivant sera randomisé dans l’autre bras. Avec 3 patients randomisés dans le même bras, il a une forte présomption. Cela peut lui permettre de différer ou annuler l’inclusion d’un patient qu’il ne souhaiterait pas voir dans le bras de randomisation certain ou probable. Cela entraînera un déséquilibre entre les groupes sur certaines caractéristiques telle que la sévérité de la maladie. Il peut aussi y avoir un mauvais respect de l’ITT ou une analyse en per protocol qui exclue plus de patients dans un groupe que dans un autre.

C’est ainsi qu’on peut interpréter des petits p comme des détecteurs à défaut de randomisation. Considérons le seuil de significativité 10^-6 et 20 caractéristiques comparées à chaque essai clinique. Si on considère le risque de grosse triche sur la randomisation à une chance sur 100 (je ne sais pas si la réalité est plutôt à 1/10, 1/100 ou 1/1000) et qu’un quart de ces triches soient suffisamment importantes sur des essais cliniques de suffisamment grande taille pour conduire à une différence significative au seuil 10^-6, alors on s’aperçoit qu’un essai sur 400 (1/400=0.0025) sera un essai avec triche de randomisation et p<10^-6 alors que 99/100×(1-(1-10^-6)^20)=1.98×10^-5 seront un essai clinique sans triche de randomisation mais avec un p<10^-6. En présence d’un petit p<10^-6 la probabilité a posteriori de triche est égale à 0.0025/(0.0025+1.98×10^-5) = 99.2%. Si on considère que la grosse triche est seulement probable à une chance sur 1000 et toujours détectable (p<10^-6) une fois sur quatre, alors la probabilité a posteriori est de 92.6%. Le doute est donc permis. Par contre, si on abaisse le seuil de significativité à 10^-9, alors même en considérant que la triche est exceptionnelle (1 chance sur 10 000) et qu’en cas de triche elle ne conduit à une différence significative au seuil 10^-9 qu’une fois sur 100, alors la proba de triche a posteriori d’un résultat significatif à 10^-9 est de 98%.

C’est pourquoi je pense que les tests statistiques restent pertinents, à condition de décaler complètement le seuil de significativité. À 10^-6 un petit p doit éveiller les soupçons et à 10^-9 il y a manifestement un problème, du moins si des tests statistiques robustes ont été appliqués. Par exemple, le test exact de Fisher fonctionne assez bien pour des petits p extrêmes (il a juste tendance à être un peu conservatif) alors que le test du chi² et le test de Student reposent sur un théorème central limite qui est d’autant plus faux qu’on s’intéresse aux queues de probabilité. Il est à noter que ces tests sont plus solides en cas de randomisation 1:1 qui reste le design le plus fréquent.

Il est à noter qu’en cas de stratification, les écarts au protocole sont bien plus faciles à voir. La plupart du temps, seul un déséquilibre minime est possible sur une variable sur laquelle une randomisation stratifiée a été appliquée. Un déséquilibre plus important que la valeur maximale théorique prouve un écart au protocole, sans le moindre doute. L’écart au protocole peut être modeste (p.e. usage d’un nombre de listes de randomisation plus élevé que précisé dans le protocole, ou usage de blocs de plus grande taille) ou pas.

Faut-il rapporter les petits p dans les tableaux 1 des essais cliniques randomisés ? Faut-il rapporter seulement les caractéristiques de tous les groupes poolés sans test ? Faut-il séparer chaque groupe sans les comparer statistiquement (comme dans le NEJM) ?

Mon opinion personnelle, c’est que l’approche du NEJM est optimale. Elle permet d’éviter les surinterprétations de petits p entre 0.01 et 0.05 de lecteurs mal avisés mais permet aussi au lecteur avisé de calculer a posteriori un petit p à 10^-9 s’il a un soupçon. Enfin, elle permet l’identification d’écarts au protocole sur la stratification de la randomisation, qui ne nécessite pas de calcul de petit p mais juste une comparaison brute de l’écart d’équilibre entre les groupes par rapport au déséquilibre théorique maximal. Enfin, un bénéfice non négligeable, à mon avis, c’est la réduction du risque de selective reporting. Si le petit p doit être calculé par les auteurs, ceux-ci sont susceptibles de supprimer la caractéristique considérée du tableau 1, par crainte que cela bloque la publication (selective reporting bias). S’ils ne testent pas, cela peut être identifié lors du reviewing (avec ou sans refus de l’article) ou lors de la lecture critique de l’article publié ; deux options meilleures que la dissimulation.

Les analyses statistiques orientées par les statistiques

Certains statisticiens orientent le choix des analyses statistiques par les résultats d’autres statistiques. Un grand classique, est illustré dans une « comparaison de deux groupes » sur une variable quantitative. Un test de normalité (p.e. Shapiro-Wilk) est d’abord fait. Si le test est significatif, alors la distribution est considérée comme non-normale et un test de Mann-Whitney est réalisé, autrement un test de Student sur séries indépendantes est réalisé. Éventuellement, un test d’égalité des variances (test de Levene ou de Bartlett) est réalisé pour décider de l’usage d’une variante de Student sans hypothèse d’égalité des variances (Aspin-Welch).

Ne faites pas ça ! Cela va à l’encontre d’un grand nombre de principes statistiques fréquentistes :

  1. La question de recherche décide de la statistique à comparer
  2. Un même protocole doit conduire à une même statistique
  3. La validité de chacun de ses tests repose sur une analyse des fluctuations d’échantillonnage qui n’est valide que lorsque la statistique est reproduite à l’identique à tous les coups

La question de recherche décide de la statistique à comparer

D’abord, le test de Mann-Whitney n’est pas un test de comparaison de moyennes. Si l’objet de l’étude est de comparer les coûts moyens de prise en charge, par exemple, on peut craindre que le test fasse une erreur de 3ème espèce, c’est-à-dire conclure à un coût moyen plus bas dans le groupe où il est plus élevé.

Le code R suivant montre un exemple dans lequel les conclusions des deux tests sont opposées :

set.seed(2020)
a=exp(2*rnorm(1000))
b=1+exp(rnorm(1000))
t.test(a,b,alternative="greater") # significativement supérieur
wilcox.test(a,b,alternative="less") # significativement inférieur

Ce phénomène paradoxal s’applique notamment à une intervention ayant un coût fixe non négligeable mais évitant des surcoûts rares et élevés. On peut, par exemple, imaginer qu’une mesure de prévention telle que l’antibioprophylaxie pour certaines chirurgies, coûte quelques euros à chaque intervention (p.e. 30 €) mais économise des milliers d’euros (p.e. 3000 € en moyenne) pour un sujet sur 50. Dans ces conditions, le coût moyen de prise en charge est plus bas avec l’intervention alors que le test de Mann-Whitney conclura à un coût plus élevé !

Cela est dû au fait que le test de Mann-Whitney ne compare pas les moyennes mais compare l’aire sous la courbe ROC à 0.50, en considérant que le test est la variable quantitative (p.e. durée de séjour) et le groupe est le diagnostic binaire. On notera au passage que le test de Mann-Whitney ne compare pas les médianes. Il m’est déjà arrivé d’observer des durées médianes de séjour identiques entre deux groupes alors que le test de Mann-Whitney montrait une différence significative.

Si on s’intéresse au coût moyen, alors on fera un test de comparaison de moyennes. Si on s’intéresse à une capacité diagnostique d’un dosage biologique, par exemple, alors on tracera une courbe ROC et le test de Mann-Whitney aura un certain sens (même si la comparaison à 0.50 n’est pas forcément très pertinente).

Un même protocole doit conduire à une même statistique

Avec la méthode de Student, on peut fournir un intervalle de confiance à 95% de la différence de moyennes. Si le protocole était répété de très nombreuses fois, dans 95% des cas, l’intervalle de confiance contiendrait la « vraie » différence de moyennes dans la population. C’est ce qu’on appelle la couverture de l’intervalle de confiance. Le test de Shapiro-Wilk étant aléatoire, de manière aléatoire, un Mann-Whitney ou un Student sera réalisé. Si la puissance du Shapiro-Wilk est de 50%, alors un Student sera fait une fois sur deux. Comment définir alors la couverture ? On peut définir la couverture conditionnelle au fait que la statistique soit générée. Auquel cas, on s’aperçoit que la statistique de Shapiro-Wilk n’est pas indépendant du résultat du test de Student, et conduit à un biais de couverture très important !

Analyse des fluctuations d’échantillonnage conditionnelles

set.seed(2020)
delta=2
v=t(sapply(1:100000, function (x) {
	a=rexp(30)
	b=delta+rexp(30)
	if (shapiro.test(c(a,b))$p.value<0.05) {
		return (c(NA,NA))
	} else {
		tt = t.test(b,a,var.equal=T)
		return(tt$conf.int)
	}
}))
mean(is.na(v[,1])) # puissance du test à 95%
# on conditionne à un Shapiro non significatif
w=v[!is.na(v[,1]),]
# On calcule les défauts de couverture
# risque borne haute 12.7% (théorie 2.5%)
mean(w[,2] < delta)
#risque borne basse 0% (théorie 2.5%)
mean(w[,1] > delta)

Ce code R montre qu’avec deux échantillons de taille 30 suivant deux lois exponentielles décalées de 2, l’intervalle de confiance à 95%, au lieu d’avoir un risque de 2.5% de surestimation et 2.5% de sous-estimation, aura un risque de 12.7% de sous-estimation et de 0% de surestimation. Cela est un biais majeur ! Si on enlève le test de Shapiro-Wilk, les risques passent à 2.4% de chaque, très proche du risque nominal de 2.5% de chaque côté. En bref, le test de Student est robuste aux écarts à la normalité, mais la séquence Shapiro->Student est complètement biaisé.

Analyse des fluctuations du petit p

Enfin, considérons les fluctuations d’échantillonnage d’un petit p issu aléatoirement d’un test de Mann-Whitney ou de Student selon le résultat d’un Shapiro-Wilk. Il est possible que les moyennes des deux groupes soient égales (première hypothèse nulle) mais que le test soit presque toujours significatif (car le test de Shapiro-Wilk est suffisamment puissant et les moyennes de rang sont inégales) conduisant à un risque alpha tendant vers 100% comme dans l’exemple ci-dessous, sous le logiciel R :

# moyennes égales...
# mais moyennes des rangs différentes
set.seed(2020)
v=sapply(1:1000, function (x) {
	a=rexp(2500) # moyenne=1
	b=rnorm(2500)+1 # moyenne=1
	if (shapiro.test(c(a,b))$p.value<0.05) {
		wilcox.test(a,b)$p.value
	}  else {
		t.test(a,b,var.equal=T)$p.value
	}
})
mean(v<0.05) # Risque alpha 99.8%

Sous la seconde hypothèse nulle, d’égalité des moyennes des rangs (hypothèse nulle de Mann-Whitney), alors il est possible d’avoir un risque alpha élevé aussi si le test de Student est souvent réalisé.

Sous l’hypothèse nulle de distribution parfaitement superposées, alors le problème est nettement moindre. Néanmoins cette hypothèse nulle peut être testée beaucoup plus efficacement par le test de Kolmogorov-Smirnov mais elle n’a aucune pertinence. Par exemple, un groupe peut simplement avoir une variance différente de l’autre, sans pour autant que la moyenne, la médiane ou la moyenne des rangs ne diffère. Un test prouvant que les deux distributions ne sont pas superposées ne permet nullement de savoir laquelle des deux distributions est la « meilleure » dès lors qu’on s’intéresse à une exposition potentiellement bénéfique (p.e. traitement) ou nocive (p.e. exposition épidémiologique).

Pour les méta-analyses

Un mélange de statistiques est presque inexploitable. Les méta-analystes sont obligés de bidouiller. Souvent ils vont approximer la médiane à la moyenne si cette première est fournie.

Au final

Le test de Student est assez robuste aux écarts à la normalité. Le test de Mann-Whitney n’est pas un test de comparaison de moyennes. La séquence des tests est une ineptie statistique. Ce problème existe aussi avec d’autres séquences de tests, comme on en retrouve, par exemple, dans les méthodes pas-à-pas, les comparaisons de modèles guidant la suite des analyses (p.e. comparaisons de plusieurs polynômes fractionnaires, tests de linéarité, comparaisons de transformations, tests « omnibus » avant de réaliser des comparaisons deux à deux, etc.)

Le choix d’une statistique doit être guidé par des considérations théoriques. Lorsqu’on s’intéresse à une différence des moyennes, plusieurs estimateurs sont possibles (Student, bootstrap) et encore une fois, ce choix doit se faire a priori plutôt qu’aléatoirement parce que chacune des procédures n’est valide que lorsqu’elle est réalisée inconditionnellement.

Modèle multivarié et essai clinique randomisé : résultat du P-hacking

Voici quelques exemples d’augmentation du risque alpha (encore appelé chance alpha) du fait de P-hacking sur les ajustements. L’analyse statistique est réalisée dans un modèle linéaire général avec petit p calculé selon une loi de Student à partir de la variance des coefficients estimée sur l’échantillon (méthode de Wald).

Le scénario consiste en une randomisation simple 1:1 de N=100 sujets (~50 par groupe) avec un outcome quantitative suivant une loi normale et 4 covariables suivant aussi des lois normales.

Le P-hacking est basé sur les libertés suivantes :

  1. Ne pas mettre la covariable du tout dans le modèle
  2. Mettre la covariable en effet quantitatif linéaire
  3. Mettre la covariable dichotomisée sur la médiane
  4. Découper la variable en tertiles, quartiles ou quintiles
  5. Et ce, séparément pour chaque variable

D’autres recodages auraient pu être envisagés, tels que le découpage à des seuils « conventionnels » (p.e. tranches de 10 ou 15 ans pour l’âge) ou l’usage de polynômes fractionnaires. Je suis resté simple. L’usage de trop de techniques n’est pas non plus habituel dans les essais cliniques randomisés.

Le P-hacking est obtenu par l’exécution de toutes les combinaisons possibles et la sélection du petit p minimal.

Dans le premier scenario, les covariables sont de réelles covariables pronostiques avec des effets beta=+1 (+1 écart-type d’outcome par écart-type de variable), beta=+0.50, beta=+0.25 et beta=+0.10 pour les 4 covariables, c’est-à-dire une covariable fortement pronostique, une moyennement pronostique et deux faiblement pronostiques. Pour un seuil de significativité bilatéral à 5%, dans ce scenario le risque (ou chance) alpha bilatéral monte à 24.7% (incertitude 24.2 à 25.2%). NB: cela correspond à une chance alpha unilatérale deux fois plus petite. C’est cette chance unilatérale qui intéresse généralement le P-hackeur.

Le second scenario est identique au premier à ceci près que le non-ajustement sur une covariable est interdit. Cela correspondait au respect d’un protocole précisant la liste des covariables d’ajustement sans en préciser le codage. La chance alpha bilatérale est à 20.9% (incertitude 20.5 à 21.4%).

Dans le troisième scenario, le non-ajustement sur une covariable est toujours interdit, mais en plus aucune des covariables n’est corrélée à l’outcome (les 4 effets sont nuls). La chance alpha redescend à 8.1% (incertitude 7.8 à 8.4%).

Le quatrième scenario est identique au troisième à ceci près que 200 sujets sont randomisés plutôt que 100. La chance alpha descend à 7.1% (incertitude 6.8 à 7.4%)

Le cinquième scenario est identique au second à ceci près que la taille d’échantillon est de 400 plutôt que 100. La chance alpha passe de 20.9% (pour N=100) à 19.8% (pour N=400, incertitude 19.4 à 20.2%).

Le sixième scenario est assez différent. Il est basé sur 20 covariables non corrélées à l’outcome. Ainsi, on compense l’absence de facteur pronostique par le nombre de variables testées. L’algorithme de P-hacking est très allégé car le nombre de combinaisons possibles est trop grand. On se contente alors de choisir en analyse trivariée (outcome ~ traitement+covariable) pour chaque covariable, le codage optimal (améliorant le petit p dans le sens de la supériorité du traitement innovant), y compris la suppression de la covariable. Une fois chacune des covariables analysée, on intègre toutes ces covariables codées comme il faut, dans le modèle le « meilleur ». Ce sixième scenario est réalisé avec N=100 et la liberté de supprimer des covariables, comme vous l’aurez compris. La chance alpha bilatérale montait à 37.4% (incertitude 36.6% à 38.2%).

Ainsi, il semblerait que le P-hacking soit efficace, même lorsque les covariables d’ajustement sont précisées dans le protocole (mais pas leur codage), mais son efficacité dépend beaucoup du degré de corrélation entre les variables pronostiques et l’outcome. Si elles sont fortement corrélées à l’outcome, alors les corrélations négatives qui peuvent aléatoirement apparaître entre le traitement innovant et le facteur pronostique vont booster le petit p d’une manière dépendante de la force pronostique. Cette force pronostique est la somme d’une partie fixe (pour les vrais facteurs pronostiques) et d’une partie aléatoire d’autant plus grande que l’échantillon est petit. C’est pourquoi, si on veut maximiser sa chance alpha, il faut rechercher de vrais facteurs pronostiques lorsque l’échantillon est grand. Pour les petits échantillons, on peut se concentrer sur de multiples variables peu ou pas corrélées à l’outcome.

Je n’ai pas analysé la chance alpha induite par l’imputation des données manquantes. Cela n’est pas évident de sélectionner des scenarii plausibles.

Il serait possible de réaliser ce travail de simulations avec une étude réelle, de telle sorte que les facteurs pronostiques aient des corrélations « de la vraie vie » avec l’outcome.

Commentaire sur « Hydroxychloroquine and azithromycin as a treatment of COVID-19: results of an open-label non-randomized clinical trial »

Billet sur un sujet d’actualité : le virus SARS-CoV-2 responsable de l’épidémie de COVID-19. Le Pr Didier Raoult, virologue français, fait la promotion du traitement par hydroxychloroquine pour le COVID-19. Notre ministre de la santé a annoncé qu’un essai clinique plus vaste soit initié, ce qui me paraît raisonnable. Que peut-on dire de la qualité méthodologique et de la fiabilité des conclusions de l’essai clinique non randomisé Marseillais ?

https://doi.org/10.1016/j.ijantimicag.2020.105949

Si je dois faire bref : il s’agit d’un essai clinique de petite taille, de méthodologie très médiocre, plutôt mal conduit, comportant de nombreux biais et dont les conclusions doivent être pris avec des pincettes. La confirmation dans un essai clinique randomisé de grande taille, avec un critère de jugement clinique, me paraît indispensable.

Critère de jugement

Le critère de jugement, purement biologique (négativation de la PCR virale) est pertinent pour un essai clinique de petite taille mais sa corrélation à un meilleur résultat clinique est loin d’être sûre, d’autant qu’il y a trois passages en réanimation, un décès et une interruption de traitement pour effet indésirable sur 26 patients traités par hydroxychloroquine, soit 5/26 = 19% (IC95% : 6.6% à 39.4%) de résultats cliniquement non satisfaisants.

Il est à noter que le critère de jugement principal décrit dans l’article est la négativation de la PCR virale à J6 alors que le protocole EudraCT parle de détection virale à J1, J4, J7 et J14, sans hiérarchiser. En bref, cela suppose des tests multiples et une correction de multiplicité des tests, et en aucun cas une détection virale à J6.

Critères d’inclusion et biais de sélection

L’essai clinique n’est pas randomisé. Le groupe intervention est constitué de patients de Marseille ayant accepté la participation à l’essai clinique à Marseille et respectant les critères d’inclusion : patient hospitalisé d’âge supérieur ou égal à 12 ans et PCR SARS-CoV-2 sur échantillon naso-pharyngé positif à l’admission, quel que soit l’état clinique. Les femmes enceintes et les sujets ayant des contre-indications à l’hydroxychloroquine étaient exclus.

Jusque là, des critères plutôt classiques, assez larges mais pragmatiques (population d’inclusion proche de la population susceptible de bénéficer du traitement), permettant une inclusion rapide de patients.

L’échantillon contrôle est beaucoup plus douteux. C’est un mélange de patients ayant de Marseille refusé l’hydroxychloroquine (nombre non précisé) et de patients recrutés à Nice, Briançon et Avignon; centres dans lesquels aucun patient ne bénéficiait du traitement innovant. Il semble étonnant que le traitement ne soit pas proposé dans les autres centres et on peut craindre que les patients des autres centres aient été recrutés de manière assez chaotique, en dehors du cadre de la recherche. Que dit le protocole posté sur EudraCT. La section E.8 est assez claire : l’essai n’est pas contrôlé (E.8.1) et il est monocentrique (E.8.3). Le protocole prévoyait 25 sujets, ce qui est cohérent avec les 26 recrutés dans le bras expérimental de l’étude (20 analysés).

Les doutes quant à la qualité du groupe contrôle sont confirmés par les données brutes fournies en annexe. Elles montrent 2 données manquantes sur le dosage de la charge virale à l’admission (J0) dans le groupe contrôle. Ces 2 sujets bénéficient d’une première PCR à J2 et à J3 respectivement. Comment est-il possible de ne pas rechercher le virus à J0 alors que ça fait partie des critères d’inclusion ? On peut craindre que les patients aient été inclus dans l’étude a posteriori, après qu’une charge virale ait été mesurée à J2 ou J3 alors que le dosage n’avait pas été fait initialement. Ce problème aurait été évité si l’étude avait bénéficié d’un e-CRF centralisé avec inclusion des patients au fur et à mesure. On peut craindre un biais de sélection secondaire important du groupe contrôle. Durant le temps où Marseille est arrivé à inclure 26 patients, les trois autres centres sont arrivés, au maximum à recruter 16 patients, sans compter la partie recrutée à Marseille. On peut donc craindre que plus de patients que ça étaient éligibles mais que seuls 16 ont été sélectionnés selon des critères que nous ne connaissons pas.

Les données publiées permettent aussi de constater que deux patients du groupe contrôle ont 10 ans et donc, ne respectent pas les critères d’inclusion (âge >= 12 ans). On constate aussi sur ces données que le statut virologique n’est pas renseigné de la même manière chez tous les patients. Pour 6 patients contrôles sur 16, le nombre de cycles (CT) de PCR est précisé lorsqu’il est en dessous de 35 alors que pour les 10 autres, la donnée est binaire « positif » ou « négatif ». Sur les 20 patients du groupe hydroxychloroquine, tous ont des mesures quantitatives. Ce défaut de standardisation des données suggèrerait l’usage d’un tableur Excel (ou LibreOffice) plutôt qu’un e-CRF pour la gestion des données, avec le défaut de traçabilité et le risque de biais de sélection secondaire qu’on peut imaginer.

Exclusions

Sur les 26 patients ayant bénéficié de l’hydroxychloroquine, 6 ont été exclus des analyses. Cela a déjà été mentionné plus haut, mais exclure un sujet décédé, trois passages en réanimation, une interruption de traitement pour effet indésirable, pour évaluer la « guérison virale » est méthodologiquement inacceptable. Le seul patient pour lequel la question se pose, c’est le sujet sorti d’hospitalisation qui avait probablement un bon état clinique. Une analyse plus pertinente aurait été une analyse en intention de traiter, en considérant tout passage en réa ou décès comme un échec. Les autres patients, notamment l’interruption thérapeutique auraient dus être suivi virologiquement comme les autres et jugés comme les autres.

Il est à noter que le résultat de l’analyse principale reste significative même en considérant que tous ces patients sont des échecs (équivalent à une positivité virale à J6). Ainsi, la mauvaise méthodologie ne remet pas forcément en cause complètement les résultats.

Données manquantes

Il y a BEAUCOUP de données manquantes et leur fréquence diffère beaucoup selon le groupe : 7/140 (5%) dans le groupe hydroxychloroquine contre 43/112 (38.4%) dans le groupe contrôle, ce qui confirme encore la non-comparabilité des groupes. Par ailleurs, ces données manquantes sont aussi présentes au moment de l’évaluation du critère de jugement principal à J6 : 1/20 (5%) dans le groupe hydroxychloroquine vs 5/16 (31.2%) dans le groupe contrôle. Une imputation LOCF a été effectuée. Même si cette imputation n’est pas aberrante, cela affecte de manière non négligeable la précision statistique (surestimée avec l’imputation simple). Par ailleurs, cela désavantage le groupe avec le plus de données manquantes (groupe contrôle) car moins de temps est laissé au patient pour se négativer. Par exemple le patient 15 (groupe contrôle) a été inclus sur un dosage virologique à J3 (écart au protocole) et a seulement eu le droit à un second dosage à J5. Par principe, il ne lui a été donné que deux jours pour se négativer là où les sujets du groupe témoin avaient généralement 6 jours.

On remarquera aussi que les détection virales « clignotent ». Certains sujets positifs, se négativent à un moment puis se re-positivent plus tard. Ainsi, 4/16 sujets clignotent dans le groupe contrôle et 4/20 sujets clignotent dans le groupe hydroxychloroquine. Cela fait douter de la qualité des prélèvements ou de la fiabilité de la PCR. Cela pose aussi le problème de la fiabilité de l’imputation LOCF puisque le sujet du groupe hydroxychloroquine, négatif aux dernières nouvelles, pour lequel on n’a pas recherché le virus à J6 est avantagé par cette imputation (étant donné qu’il aurait pu se re-positiver) tandis que les 5 sujets du groupe contrôle ont été désavantagés par la non-mesure de la charge virale à J6 puisqu’ils étaient positifs aux dernières nouvelles.

On peut encore remarquer que le groupe contrôle se divise en deux : les sujets avec quantification (n=6) et les sujets sans quantification (n=10). Il est possible que ces sujets avec quantification corresponde aux sujets de Marseille (ce ne serait pas étonnant, vu qu’il y a 4 enfants asymptomatiques sur six patients, pour lesquels on peut comprendre que les parents aient refusé la participation à un essai clinique sur un traitement innovant alors que tout allait bien pour leur petit bout de chou) ou alors il s’agit au moins de sujets suivis avec une « quantification précise » de la PCR. Considérons donc le groupe contrôle. Sur les 6 sujets avec quantification précise, 5 ont au moins une virologie négative à un moment où à un autre contre 0 sur les 10 sujets sans quantification précise (p=0.0014 selon un test exact de Fisher). Comme il s’agit d’une analyse post hoc, ce résultat est à prendre avec des pincettes (au même sens que l’analyse très douteuse sur l’Azithromycine). Une partie de la différence est explicable par le nombre de données manquantes extrêmement grand chez les sujets sans suivi rigoureux (38/70=54%) baissant les chances d’avoir un prélèvement négatif dans le long. Néanmoins, on peut douter de la comparabilité des prélèvements des patients sans quantification précise : le nombre de cycles de PCR était-il le même ? Les conditions de prélèvement étaient-elles les mêmes ?

Réanalyse avec nouvelle gestion des données manquantes

D’abord, les six sujets du groupe hydroxychloroquine exclus ont été réintroduits avec l’imputation suivante : le décédé et les sujets passés en réa ont été considérés comme d’évolution favorable (positif à J6). Les deux autres sujets (effet indésirable et sortie d’hospitalisation) ont été analysés sur la base des données disponibles, les prélèvements ultérieurs étant considérés comme « non faits ».

La première ré-analyse a été faite avec une imputation multiple par équations chaînées selon la méthode du package mice. Les variables utilisées pour l’imputation étaient : le groupe de traitement et chacun des dosages binaires (positif ou négatif) de D0 à D6. Chaque variable était imputée dans un modèle de régression logistique expliqué par toutes les autres variables. Un total de 500 jeux d’imputation a été réalisé. Une régression logistique bivariée expliquant la positivité à J6 selon le groupe de traitmeent était estimé par le maximum de vraisemblance pour chaque jeu d’imputation. Les log-odds-ratio ont été poolés avec la méthode de Rubin, c’est-à-dire, globalement en additionnant la variance intra à la variance inter. La méthode de Wald, prenant en compte le nombre de degrés de liberté a été finalement utilisé pour estimer le degré de significativité de la statistique poolée ainsi que son intervalle de confiance.

La fraction d’information manquante (FMI) était estimée à 0.26. L’odds ratio de l’effet du traitement était estimé à 0.143 (IC95% : 0.023 à 0.884, p=0.037). Cette estimation est très fragile car les conditions de validité asymptotiques des méthodes utilisées ne sont pas atteintes.

Une autre méthode de calcul, c’est simplement un test exact de Fisher en excluant les sujets avec une évaluation manquante à J6. On obtient alors un odds ratio à 0.18 (IC95% : 0.016 à 1.15, p=0.064).

Il semble donc que la tendance statistique persiste avec ces analyses de sensibilité. Néanmoins, la fiabilité de résultats en présence d’une mauvaise qualité de données reste toujours sujette à caution.

Deuxième ré-analyse : comparaison à des données historiques

Problème de multiplicité des tests

La comparaison entre les groupes est répétée tous les jours de J1 à J6, soit six tests pour la comparaison hydroxychloroquine vs contrôle. Cela peut paraître beaucoup mais ne me préoccupe pas tant que ça. D’abord, l’analyse principale est annoncée comme la comparaison à J6. C’est néanmoins douteux car ça ne correspond pas au protocole EudraCT. On peut alors craindre un choix a posteriori du critère de jugement principal et donc craindre du P-hacking. Sans préjuger d’autres formes de P-hacking, la sémiologie du P-hacking par multiplicité des tests n’est pas retrouvée. Cette sémiologie correspond à un résultat à la limite de la significativité statistique (p entre 0.01 et 0.05), sans cohérence globale telle qu’une différence qui aurait un p=0.03 à J4 mais un p=0.30 à J3 et un p=0.40 à J5. Dans cette étude, les différences entre les deux groupes sont retrouvées de manière cohérente avec un écart progressivement croissant et une forte significativité. Le bémol, c’est que les écarts se resserrent dans les analyses de sensibilité que j’ai réalisées. C’est donc moins robuste qu’il paraît.

Par contre, l’analyse en sous-groupe sur l’Azithromycine+Hydroxychloroquine est très douteuse. Elle n’était absolument pas planifiée dans le protocole EudraCT et pourtant se retrouve dans le titre de l’article, correspond à un sous-groupe minuscule (n=6). Déjà, le test statistique employé est incorrect. C’est un test global sur les trois groupes, répondant à la question : y a-t-il un au moins des trois groupes dont le pourcentage de négativation à J6 diffère des autres ? Pour conclure au bénéfice de l’Azithromycine en addition à l’Hydroxychloroquine et en considérant que cette dernière a déjà fait preuve de son efficacité (avec tous les bémols mentionnés au dessus), alors il faudrait comparer le groupe Hydroxychloroquine seul à Azithromycine+Hydroxychloroquine. Avec le test exact de Fisher, le taux de succès à J6 de 8/14 n’est pas significativement différent de 6/6 (p=0.115). Par ailleurs, d’autres analyses en sous-groupes étaient possibles : par exemple Azithromycine seule vs contrôle ou Azithromycine±Hydroxychloroquine vs Autre. Il aurait pu aussi y avoir d’autres analyses en sous-groupes sur la prise en charge… On ne peut pas savoir combien ont été faites étant donné qu’aucun protocole publié ne fournit d’information dessus.

D’une manière générale, on évite les analyses post hoc en sous-groupes dans les essais cliniques, mais c’est encore plus vrai lorsqu’ils sont de toute petite taille, car la puissance est proche du risque alpha, conduisant à un risque de fausse découverte très élevé.

Ajustements statistiques ?

Les effectifs sont tellement petits qu’il est difficile de faire une comparaison statistique des caractéristiques des patients. Ces comparaisons ont tendance à être non significatives alors que des différences majeures existent entre les groupes. Ces différences peuvent être expliquées par les fluctuations d’échantillonnage, ou être expliquées par une sélection différente des patients, mais faire la part des choses est très difficile. On a quand même quelques indices pour penser que la population traitée a une sévérité de la maladie un plus élevée : trois passages en réanimation et un décès sur 26 patients (versus zéro sur 16 ?). Une moyenne d’âge plus élevée (52.1±18.7 vs 37.3±24 ans, p=0.06) chez les sujets non exclus. La différence aurait peut-être été encore plus forte si on avait inclus les patients passés en réa et décédés puisque ces risques touchent quand même nettement plus les personnes âgées. Archétype de la variable impossible à comparer : la proportion d’infections respiratoires basses (LRTI) est de 30% dans le groupe hydroxychloroquine vs 12.5% dans le groupe contrôle, mais le nombre de sujets est tellement bas que ça peut être complètement expliqué par le hasard (p=0.30).

Certains pourraient proposer des ajustements statistiques. C’est difficile avec des effectifs aussi faibles. On peut juste se demander dans quel sens va le biais. La charge virale semble plutôt positivement corrélée à la sévérité (https://doi.org/10.1016/S1473-3099(20)30232-2) et la clairance semble être plus lente dans les formes sévères. Ainsi, le groupe hydroxychloroquine est a priori désavantagé par ce biais de sélection différentiel des formes plutôt sévères. Cela ne remet donc pas en cause les résultats.

Nombre de sujets nécessaires

L’article comporte un calcul de nombre de sujets nécessaires comme si l’étude avait été planifiée avec un groupe contrôle dès le départ. Cela semble être un calcul a posteriori. Je n’ai pas de preuve absolue mais le fait que le protocole EudraCT soit clairement conçu comme un essai à un seul bras et est un premier indice. Le second indice est le fait que le groupe contrôle soit construit d’une manière très étrange, mélangeant des refus de participation à l’étude (j’espère qu’ils ont quand même signé un consentement pour avoir des PCR tous les jours, sinon le CPP pourrait ne pas approuver) et des patients des centres non Marseillais pour lesquels le traitement n’était pas une option. L’inclusion de ce dernier groupe de patients non Marseillais paraîtrait éthiquement douteux si les patients avaient des prélèvements sanguins de PCR spécifiquement pour la recherche sans avoir la possibilité de bénéficier du traitement innovant ; mais il paraît possible qu’ils aient simplement été inclus comme une série de cas rétrospective à partir des dossiers médicaux et qu’en conséquence ils n’aient eu aucun examen ou traitement en dehors de la routine.

Au total

Le biais le plus préoccupant de cette étude, à mon avis, provient du suivi et surtout de la sélection des contrôles en dehors de Marseille. Certains de ces sujets n’ont pas eu de PCR à baseline mais ont quand même été inclus. De très nombreuses données manquent sur leur PCR virales. Leur PCR est considérée comme « positive » ou « négative » sans quantification, pour la plupart. On peut se permettre de douter du fait que les prélèvements ont été réalisés, cultivés et analysés de la même manière qu’à Marseille. Mais le pire, c’est qu’on ne sait pas quels sujets étaient éligibles dans le groupe contrôle mais n’ont pas été inclus dans l’étude.

Sans ce problème de groupe contrôle, le résultat principal à J6 aurait une certaine robustesse malgré les données manquantes et la multiplicité potentielle des tests (puisque le CJP n’est pas celui du protocole). Du fait de ce groupe contrôle douteux, j’émets des réserves quand à la reproductibilité de l’analyse principale.

Cela ne veut pas dire que le traitement est inefficace ou que l’article n’apporte rien. Le résultat est très fragile mais fournit quand même une piste intéressante, à mon avis. Je me range du côté de la communauté scientifique qui propose de réaliser un essai clinique randomisé de grande ampleur afin de statuer sur l’efficacité clinique de l’hydroxychloroquine avec une méthodologie robuste.

Ajout post hoc : comparaison à des données historiques

Comme décrit avant, la fragilité principale de l’étude provient de la sélection du groupe contrôle sans dosage quantitatif précis (c’est-à-dire avec juste POSITIF ou NEGATIF comme résultat de la PCR). On peut craindre qu’il s’agisse de données de routine rétrospectivement collectées, avec un biais de sélection. On peut notamment craindre qu’il y ait eu l’exigence de la présence d’un prélèvement à J5 ou J6 pour inclure le patient. Cela risque d’exclure les patients ayant eu une négativation précoce pour lesquels il n’y a plus de raison de continuer les prélèvements au delà d’un ou deux résultats négatifs.

C’est pourquoi il est encore préférable de comparer l’échantillon hydroxychloroquine à un échantillon indépendant, de données historiques, dont la rigueur de recueil est meilleure. Prenons comme référence l’article Viral dynamics in mild and severe cases of COVID-19. L’article ne présente que les différences de Ct par rapport au prélèvement initial (delta(Ct)). On peut néanmoins comprendre, de la figure, que les sujets n’ont plus de prélèvements une fois négatifs. Dans les formes légères la majorité des sujets est négatif à J8. À J6, il y a nettement plus de positifs : environ la moitié des patients ayant une forme légères. Les formes sévères (dont 23/30=77% sont passés en unité de soins intensifs) par contre, nécessitent nettement plus de temps (> 15 jours pour la négativation). Si on considère l’échantillon des 20 patients dont les cas réanimatoires a été exclu, la majorité doit ne pas être sévère. Un taux de 50% de négativation à J6 paraît plausible. Le taux observé est de 14/20=70% (IC95% : 46 à 88%), compatible avec le taux de 50% ; mais il y a une tendance statistique à un taux plus élevé.

La conclusion reste la même. Un effet biologique paraît plausible mais il y a un doute non négligeable. Le bénéfice clinique reste très douteux. Avant d’arroser la population d’hydroxychloroquine, je conseillerai d’inclure des patients dans l’essai clinique de grande taille ou d’en attendre le résultat.