Ajustement sur effet centre

Imaginez que vous vous intéressiez à comparer des techniques médicales ou chirurgicales pour lesquels il existe à la fois une grande variance de pratique inter-centre mais aussi une variance de pratique intra-centre. Que faire de l’effet centre ? Faut-il ajuster dessus ?

Il existe deux approches orthogonales, toutes deux soumises à des biais bien différents.

  1. La comparaison intra-centre, obtenue par analyse de la corrélation individuelle entre l’exposition (traitement) et l’outcome avec un ajustement sur l’effet centre
  2. La comparaison inter-centre, écologique, par analyse de la corrélation de l’exposition moyenne du centre avec l’outcome moyen du centre. On considère alors qu’il n’y a qu’une seule observation par centre, tous les patients de chaque centre étant aggrégés.

La première approche, est soumise à un possible biais d’indication, qui est parfois majeur, le traitement n’étant pas donné au hasard mais pouvant dépendre très fortement des caractéristiques du patient et de son pronostic. La deuxième approche ne souffre pas du biais d’indication car tous les patients d’un même centre sont assignés à la même valeur synthétique (taux moyen d’exposition) peu importe leurs caractéristiques. Cependant, elle souffre de problèmes de confusion au niveau des centres. Par exemple, on peut craindre que les centres publics de grande taille accueillent une population plus fragile, présentant un plus grand nombre de comorbidités, et aient aussi une attitude thérapeutique différente. Ces facteurs de confusion sont totalement anéantis par les comparaisons intra-centres (première approche). Dans la deuxième approche, il reste possible de procéder à des ajustements sur les caractéristiques du centre (public/privé, nombre de lits, etc.) ainsi que sur les caractéristiques aggrégées des patients (niveau de comorbidité moyen des patients). Cependant, seuls les facteurs de confusion observables et mesurés sans erreur de mesure majeure sont contrôlables alors que l’approche de comparaison intra-centre fait disparaître même les différences inobservables entre centres.

D’autres aspects pragmatiques peuvent aider à décider d’une méthode d’analyse, de l’autre ou des deux.

Si la variance de pratique inter-centre est trop faible (p.e. les centres prescrivant peu la thérapeutique sont à 45% de prescription alors que les centres la prescrivant beaucoup sont à 50%), la puissance de la deuxième approche peut être trop faible, avec un risque que l’effet soit noyé dans les biais dans le pire des cas ; en effet, l’effet du traitement sera dilué sans forcément que les biais soient diminués, conduisant à un rapport effet/biais défavorable. Au contraire, si la variance inter-centre est forte alors que la variance intra-centre est faible (p.e. les centres prescrivant peu la thérapeutique sont à 5% alors que les centres la prescrivant beaucoup sont à 95%), alors que l’approche de comparaison intra-centre sera certainement soumise à un biais d’indication incontrôlable, les rares patients n’ayant pas la pratique habituelle du centre correspondant possiblement à des indications formelles ou des contre-indications formelles de la thérapeutique ; sans compter le problème de la faible puissance de la première approche.

Un autre cas de figure doit rentrer en ligne de compte : le nombre de centres. Si celui-ci est élevé, l’approche écologique est tout à fait possible. S’il n’y a que deux ou trois centres, elle peut être impossible. On ne pourra plus estimer empiriquement la variance inter-centre de l’outcome moyen, la variance inter-centre de l’exposition moyenne et la covariance entre les deux, comme cela est possible avec l’estimation d’un coefficient de corrélation de Pearson. En effet, avec n=2, on ne peut pas estimer un coefficient de corrélation de Pearson. La comparaison inter-centre sera alors faite avec une méthodologie ici-ailleurs dans laquelle la différence du résultat moyen entre les centres sera entièrement attribuée aux fluctuations d’échantillonnages intra-centre et à la différence du taux moyen d’exposition. C’est alors que des différences de caractéristiques des centres ou de leur recrutement qui ne correspondent même pas à des facteurs de confusion vont engendrer un biais dans l’estimation de l’effet du traitement. Par exemple, si, par hasard, le centre ayant un niveau élevé de la prescription s’avère avoir une population particulièrement fragile en raison du voisinage d’une clinique absorbant l’activité des patients les plus stables, alors cette différence sera attribuée à la différence d’exposition au traitement. Une telle différence conjoncturelle serait incluse dans le calcul empirique de la variance inter-centre dans une analyse écologique intégrant un grand nombre de centres. Une étude ici-ailleurs bicentrique sera alors plus sujette au biais qu’une grande étude écologique et l’approche d’ajustement sur l’effet centre pourra éventuellement devenir plus intéressante.

Selon la nature du critère de jugement, la problématique pourra différer notablement, notamment pour une étude ici-ailleurs. La durée d’hospitalisation, par exemple, est soumise à des différences inter-centres majeures explicables par des différences organisationnelles sans forcément de rapport avec les caractéristiques des patients. Une étude ici-ailleurs comparant la durée d’hospitalisation est soumise à un biais de comparabilité initiale des centres totalement inacceptable. Au contraire, une comparaison intra-centre la durée d’hospitalisation est possible, l’organisation du service étant alors superposable pour les différents patients suivis dans un même centre sur une courte période de temps.

Les différences organisationnelles entre services devraient affecter moins la mortalité à 30 jours que la durée d’hospitalisation. Cette mortalité à 30 jours dépendra certainement des caractéristiques initiales des patients, mais, dans bon nombre de situations, on peut s’attendre à ce que les différences de caractéristiques moyennes des patients entre les centres soient plus faibles que les différenes intra-centre mais inter-traitement.

Enfin, dans certaines situations où les deux approches, de comparaison intra-centre ou inter-centre paraissent toutes deux aussi biaisées, il est possible de réaliser les deux analyses ; formellement on pourra définir une analyse principale et une analyse de sensibilité, mais l’interprétation des résultats sera conjointe. Si les deux analyses fournissent des résultats cohérents, on pourra être relativement rassuré, se disant qu’il suffit qu’une des deux analyses ne soit pas biaisée de manière majeure pour que les conclusions soient fiables. Au contraire, en cas de résultat discordant, cela montrera l’existence de biais majeurs dans au moins une des deux analyses ; considérant que les deux paraissent également biaisées initialement, on ne peut alors savoir laquelle est correcte, et les résultats sont alors très douteux.

L’objectif de ce billet était de faire comprendre que l’idée selon laquelle un ajustement ne peut qu’enlever du biais et pas en ajouter, est simpliste. L’ajustement peut créer un biais plus grand qu’il n’en corrige, et il est parfois préférable de faire une forme d’anti-ajustement sur l’effet centre par aggrégation des résultats au niveau du centre, dans une approche écologique plutôt qu’individuelle. Les choix statistiques dépendent de nombreux facteurs, et notamment d’une évaluation subjective de l’ampleur des biais. En corrigeant des biais on en crée parfois d’autres. Il faut alors se demander si les biais créés sont plus grands que les biais corrigés. Il n’existe en pratique pas de quantification objective de ces biais, car généralement, un biais quantifiable est rectifiable et rectifié. Les biais restant sont alors seulement quantifiables subjectivement. La connaissance du sujet de recherche est alors nécessaire à un jugement subjectif pertinent. Cependant, cette connaissance se heurtera aux limites des connaissances et de la compréhension humaine.

En tout cas, l’ajustement ce n’est pas automatique !

Modèles à risques compétitifs

Introduction

Le modèle de Cox, et l’estimateur de Kaplan-Meier servent classiquement aux analyses de survie. Ils reposent sur la notion d’événement, tel que le décès, dont on puisse dater précisément la survenue. Ces modèles atteignent leurs limites lorsqu’on souhaite analyser plusieurs événements mutuellement exclusifs (compétitifs).

Par exemple, on peut se demander si le cancer (survenant souvent aux alentours de 60 ans) est un facteur de risque de pemphigoïde bulleuse (survenant généralement au-delà de 75 ans). On peut analyser cela sur une cohorte de patients nouvellement diagnostiqués d’un cancer (10.1111/bjd.17197) avec un groupe contrôle sans cancer apparié sur l’âge et le sexe. Deux risques compétitifs existent alors:

  1. Le risque de pemhigoïde bulleuse
  2. Le risque de décès

La survenue d’un décès survenant avant le développement de la pemphigoïde bulleuse, empêche totalement la possibilité d’observer le développement de la pemphigoïde bulleuse chez ce patient.

Une manière de théoriser ce problème, c’est de considérer qu’il existe deux courbes latentes (c’est-à-dire non observables directement) de risque, correspondant aux deux événements (décès et pemphigoïde bulleuse). Plutôt que d’observer chacune des deux courbes, on observe seulement le premier des deux événements qui survient, chez chacun des sujets. On souhaite alors reconstruire les deux courbes de risques latentes (décès et pemphigoïde bulleuse) à partir de la courbe combinée observée.

La non-identifiabilité de ces courbes a été démontré par Tsiatis en 1975 (doi: 10.1073/pnas.72.1.20). En effet, plusieurs combinaisons possibles de courbes de risques latentes peuvent conduire à la même courbe combinée, dont notamment, une paire de courbes dans lesquelles les risques sont totalement indépendants. Le problème est donc insoluble.

Cela n’a pas empêcher deux solutions d’émerger : l’approche de Kalbfleisch & Prentice (doi: 10.2307/2530374) en 1978 et l’approche de Fine & Gray en 1999 (doi: 10.1080/01621459.1999.10474144).

Quelles sont les forces et les limites de chacune des deux approches ? Quand utiliser l’une et quand utiliser l’autre ? C’est à ces questions que nous allons tenter de répondre.

Faux problème : analyse de la survie

La théorie des risques compétitifs nous dit que la survenue du premier des événements empêche la survenue du second évenement. Ainsi, la survenue d’un décès empêche l’observation du développement de la pemphigoïde bulleuse. Cependant, dans l’autre sens, la théorie des risques compétitifs nous dit qu’il serait impossible d’observer le décès après la pemphigoïde bulleuse. Cela est ridicule, et il est totalement inutile, et même contre-productif de recourir aux modèles à risques compétitifs pour analyser le risque de décès dans une population susceptible de développer une pemphigoïde bulleuse. Il suffit d’ignorer complètement les événements de pemphigoïde bulleuse et d’utiliser les modèles classiques d’analyse de survie pour l’événement du décès : modèle de Cox, estimateur de Kaplan-Meier, etc.

Le vrai problème, c’est l’analyse du risque de développer la pemphigoïde bulleuse dans une population dans laquelle des sujets peuvent décéder avant d’avoir développé la pemphigoïde bulleuse.

Approche de Kalbfleisch & Prentice

L’approche de Kalbfleisch & Prentice (doi: 10.2307/2530374) est aussi connue sous le nom de « cause-specific hazards » et correspond à l’usage d’une fonction de risque instantané (hazard function) au temps t, conditionnelle à la survie sans le moindre événement jusqu’à ce temps t. Pour ceux qui voudraient aller plus loin, regardez la deuxième formule, page 543 de Biometrics, Vol. 34, No. 4 (Dec., 1978).

En reprenant l’exemple du décès et de la pemphigoïde bulleuse, l’inférence de ce modèle est la même que celle que l’on peut obtenir à partir des modèles usuels (modèle de Cox, estimateur de Kaplan-Meier, modèle paramétrique de Weibull) en censurant sur le décès pour analyser la pemphigoïde bulleuse. On n’a donc pas besoin de procédure spéciale d’un logiciel statistique pour appliquer cette approche. Les procédures usuelles font l’affaire. Il suffit de remplacer les décès par des censures à la date du décès.

Toujours dans le contexte de la pemphigoïde bulleuse avec risque compétitif de décès, deux interprétations de l’approche de Kalbfleisch & Prentice sont possibles :

  1. On modélise le risque de la pemphigoïde bulleuse chez les sujets n’étant pas décédés avant le moment où le risque de pemphigoïde est analysé. C’est un risque conditionnel à la survie. On comprend alors que le biais du survivant s’applique, et qu’une exposition qui éliminerait préférentiellement des sujets prédisposés à la pemphigoïde bulleuse paraîtrait être un facteur protecteur de pemphigoïde bulleuse.
  2. On modélise la courbe de risque latente de la pemphigoïde bulleuse sous l’hypothèse d’indépendance du risque de décès et du risque de pemphigoïde bulleuse.

L’interprétation N°2 étant fondée sur une hypothèse invérifiable et généralement peu plausible, je conseille très vivement l’interprétation N°1 : risque de pemphigoïde bulleuse conditionnel à la survie.

Cette interprétation N°1 est sujette au biais du survivant, ce qui peut être préoccupant pour l’interprétation d’un lien causal entre cancer et le développement d’une pemphigoïde bulleuse. C’est beaucoup moins préoccupant pour une approche pronostique. Si, par exemple, un patient ayant survécu à un cancer développé à 60 ans, demande, à l’âge de 75 ans, s’il a à craindre la pemphigoïde bulleuse, il serait pertinent, dans la réponse, de prendre en compte cet antécédent, y compris sur son possible effet d’élimination des sujets les plus fragiles à la pemphigoïde bulleuse. Le patient est un survivant du cancer, il faut donc lui fournir des statistiques concernant les survivants ! En pratique, le problème du pronostic se pose rarement comme ça, mais ce n’est qu’un exemple pour illustrer la divergence d’interprétation selon que l’on s’intéresse à la causalité ou au pronostic.

Il est à noter qu’avec l’approche de Kalbfleisch & Prentice, on peut calculer les chances d’être en vie sans pemphigoïde bulleuse au temps T, en multipliant les chances de survie brute (en ignorant le risque de pemphigoïde bulleuse) au temps T par la chance de ne pas faire de pemphigoïde bulleuse chez les sujets ayant survécu au temps T. Si on note S(t) la probabilité de survie sans décès au temps t, et SPB(t) la probabilité d’avoir fait une pemphigoïde bulleuse à un temps inférieur ou égal à t chez les survivants au temps t, selon Kalbfleisch & Prentice, alors on peut calculer trois probabilités dont la somme sera égale à 1 :

  1. Probabilité d’être décédé au temps t, après avoir fait une pemphigoïde bulleuse ou pas : 1-S(t)
  2. Probabilité d’être en vie sans pemphigoïde bulleuse au temps t : S(t)*SPB(t)
  3. Probabilité d’être en vie avec pemphigoïde bulleuse au temps t : S(t)*(1-SPB(t))

Approche de Fine & Gray

L’approche de Fine & Gray (doi: 10.1080/01621459.1999.10474144) est celle des « subdistribution hazards ». Si on reprend l’exemple de la pemphigoïde bulleuse et du décès, le concept est basé sur le fait que l’on peut définir trois probabilités dont la somme fait 1, à tout temps t :

  1. Probabilité d’être survivant sans pemphigoïde bulleuse au temps t
  2. Probabilité d’être décédé au temps t, avant d’avoir eu le temps de faire une pemphigoïde bulleuse
  3. Probabilité d’avoir fait une pemphigoïde bulleuse, avec ou sans décès ultérieur, au temps t

Avec ces trois probabilités, on peut tracer trois courbes qui, cumulées, font 100% à tous les temps t.

Les sujets de type N°1 n’ont fait aucun des deux événements au temps t. Les sujets du type N°2 ont fait l’événement de décès en premier, et donc, ne peuvent pas avoir fait de pemphigoïde bulleuse. Les sujets de type N°3 ont fait l’événement de pemphigoïde bulleuse, et donc, ne peuvent pas décéder selon la théorie. En réalité, ils peuvent décéder, mais c’est complètement ignoré par le modèle, parce que la théorie des risques compétitifs, nous dit que dès qu’un des deux événements est survenu (décès ou pemphigoïde bulleuse), l’autre événement ne peut plus être observé (respectivement, pemphigoïde bulleuse ou décès). Pour la théorie des subdistributions hazards, ça va plus loin. On considère que dès qu’un événement survient (décès ou pemphigoïde bulleuse ), l’autre événement ne peut plus survenir ( respectivement, pemphigoïde bulleuse ou décès) et ne doit alors plus jamais être comptabilisé dans les subdistribution hazards.

Pour mieux comprendre cela, il faut voir comment on peut estimer ce modèle. Dans le cas simple où tous les patients sont suivis jusqu’au décès ou à la pemphigoïde bulleuse, c’est-à-dire, qu’il n’y a aucune censure administrative et aucun perdu de vue, alors le modèle de Fine & Gray décrivant le risque de pemphigoïde bulleuse est équivalent à un modèle de Cox dans lequel on recoderait les sujets décédés comme étant indemnes de la pemphigoïde bulleuse jusqu’à un temps virtuellement infini, en les censurant à une valeur supérieure à tous les délais qui ont été observés dans la base de données. Par exemple, si le suivi maximal est de 30 ans, on peut considérer que les sujets décédés sont survivants sans pemphigoïde bulleuse à 1000 ans en les censurant à 1000 ans. C’est un tout petit peu plus compliqué s’il y a de vraies censures, à cause de perdus de vue ou sujets partiellement suivis, mais l’idée restera la même : les sujets qui sont décédés deviennent invulnérables à la pemphigoïde bulleuse mais continuent de compter comme sujets « à risque de pemphigoïde bulleuse » dans les courbes de survie.

Quelles sont les conséquences pratiques de tout ça ? Comment peut-on interpréter les résultats ?

La première chose, c’est qu’il ne faut pas du tout utiliser ce modèle pour comparer le risque de décès entre les groupes, car le modèle considère que tout patient qui a fait une pemphigoïde bulleuse est immortel au-delà de sa pemphigoïde bulleuse. De manière mécanique, toute exposition qui augmente le risque de pemphigoïde bulleuse aura l’air de diminuer le risque de décès car en déclenchant des pemphigoïdes bulleuses, il déclence des immortalités. Ce n’est pas grave, parce que l’analyse du risque de décès ne nécessite pas de modèle à risque compétitif comme précisé dans la section « Faux problème : analyse de la survie ».

La deuxième chose, c’est la manière d’interpréter le risque de pemphigoïde bulleuse dans ce modèle. Tout sujet décédé rentrant automatiquement en catégorie N°2 (décédé avant pemphigoïde bulleuse), il sera considéré comme invulnérable au risque de pemphigoïde bulleuse et ne pourra plus du tout rentrer dans la catégorie N°3 (ayant fait la pemphigoïde bulleuse avant de décéder). En conséquence, toute maladie augmentant de manière importante le risque de décès, comme le cancer, aura l’air de protéger de manière majeure du risque de pemphigoïde bulleuse. Par exemple, si le cancer tue 50% des patients rapidement, et que le risque de pemphigoïde bulleuse est le même chez les survivants de cancers que chez les sujets qui n’ont jamais fait de cancer, le cancer aura l’air d’être un facteur protecteur divisant par deux le risque de pemphigoïde bulleuse. D’une certaine manière, ce n’est pas faux. L’euthanasie est l’ultime méthode de prévention de toutes les maladies.

Ce modèle a du sens pour répondre à la question que peut poser le patient atteint d’un cancer : « docteur, est-ce que vous croyez que j’ai des chances de voir mes petits enfants ? ». En effet, même si le décès provoqué par le cancer, ne va pas empêcher les petits enfants de naître, il va empêcher le patient de voir ses petits enfants…

Quelle approche utiliser dans la pratique ?

Avec la pemphigoïde bulleuse, il est évident que le modèle de Fine & Gray est très trompeur. Il donne l’illusion d’un effet protecteur majeur du cancer (Hazard Ratio = 0.47) sur la pemphigoïde bulleuse (10.1111/bjd.17197) parce que celui-ci tue les patients avant qu’ils aient le temps de faire une pemphigoïde bulleuse. Le modèle de Kalbfleisch & Prentice est beaucoup plus pertinent, même s’il va être soumis à un léger biais du survivant. Si le cancer tue de manière indistincte les sujets à faible et à fort risque de pemphigoïde bulleuse et ne provoque ni ne protège de la pemphigoïde bulleuse, alors, le hazard ratio de Kalbfleisch & Prentice sera égal à 1, alors que le modèle de Fine & Gray aura un Hazard Ratio égal à 0,50 si le cancer tue un patient sur deux rapidement, avant qu’il ait le temps de développer une pemphigoïde bulleuse.

Néanmoins, il existe des situations ou le modèle de Fine & Gray est bien préférable ! Lorsque les deux événements sont de nature complètement opposée. Par exemple, on peut souhaiter analyser le délai avant guérison d’une maladie. Que faire des patients décédés ? Le modèle de Fine & Gray va considérer que les sujets décédés ne vont jamais guérir. C’est certainement une hypothèse plus raisonnable que le Kalbfleisch & Prentice, qui, selon le point de vue, va considérer que les chances de guérison après décès sont les mêmes que chez les survivants, ou s’intéresser seulement au risque de guérison chez les survivants.

Voilà donc des étapes simple pour choisir le modèle.

  1. Catégoriser les événements en bénfiques (p.e. guérison, sortie d’hospitalisation, amélioration d’un état clinique) et nocifs (p.e. apparition ou aggravation d’une maladie, décès)
  2. Si les deux événements analysés sont tous deux nocifs, alors utiliser l’approche de Kalbfleisch & Prentice
  3. Si l’un des événements est bénéfique (p.e. guérison) alors que l’autre est nocif (p.e. décès), utiliser l’approche de Fine & Gray, dans laquelle on considèrera que le sujet ayant eu l’effet nocif n’aura jamais l’effet bénéfique et le sujet ayant eu l’effet bénéfique n’aura pas l’effet nocif.

Personnellement, je n’ai jamais rencontré de situation où les deux événements étaient bénéfiques, c’est pourquoi je n’ai pas listé cette possibilité dans le choix. Le risque compétitif est habituellement le décès.

Par ailleurs, même si on utilise un modèle à risques compétitifs, il faut seulement l’utiliser pour l’événement qui est caché par l’autre. La pemphigoïde bulleuse, par exemple, est cachée par le décès, mais dans l’autre sens, la pemphigoïde bulleuse ne cachera pas un décès qui peut toujours être observé ultérieurement (cf section « faux problème : analyse de la survie »).

Quelques exemples d’événements opposés

Comme précédemment décrit, le Fine & Gray est adapté aux événements opposés. Le premier cas, c’est l’analyse du délai avant guérison d’une maladie dans laquelle on peut considérer que les sujets décédés ne guériront jamais. Dans l’analyse de la survie lors d’une hospitalisation, sur une cohorte de patients pour lesquels il n’existe aucun suivi après la sortie d’hospitalisation, comme on peut parfois rencontrer dans des études rétrospectives sur données hospitalières, on peut considérer que les sujets qui sortent vivants de l’hôpital ne décèderont jamais (Fine & Gray)… Évidemment, certains décèderont, mais si on fixe un horizon temporel relativement limité (p.e. 30 à 90 jours), c’est une bonne approximation ; certainement bien meilleure que de considérer que les sujets qui sortent d’hôpital ont le même risque de décéder à court terme que les sujets encore hospitalisés au même moment (Kalbfleisch & Prentice). C’est dans ce cas de figure qu’un article du BMJ (doi : 10.1136/bmj.m1966) a utilisé à juste titre le modèle de Fine & Gray. Personnellement, j’aurais fait un peu différemment ; plutôt que de calculer des hazard ratio de décès intra-hospitaliers avec sortie comme risque compétitif, j’aurais calculé des risques relatifs ou différences absolues de risque de décès à J90, comptabilisant les décès intra-hospitaliers, mais aussi les décès durant de possibles réhospitalisations. Cependant, leur méthode n’est pas incorrecte.

Situation avec trois événements ou plus

Comment peut-on analyser trois risques à la fois, comme le décès, le risque d’accident vasculaire cérébral (AVC) et le risque d’infarctus du myocarde (IDM) ?

La théorie des risques compétitifs selon laquelle l’observation d’un événement empêche l’observation des autres est fausse dans ce cas là. On peut tout à fait observer un AVC après un IDM ou vice-versa. Il n’est alors pas nécessaire de modéliser les trois événements (AVC, IDM et décès) dans un seul modèle. Je conseillerai d’analyser le risque d’AVC avec le décès comme risque compétitif avec l’approche de Kalbfleisch & Prentice, c’est-à-dire, en censurant sur le décès. Un deuxième modèle servirait à analyser le risque d’IDM avec le décès comme risque compétitif selon Kalbfleisch & Prentice.

D’un point de vue théorique, Kalbfleisch & Prentice et Fine & Gray peuvent tout deux supporter les situations où il existe au moins trois risques vraiment compétitifs. Avec Kalbfleisch & Prentice, pour analyser un événement, on censure sur tous les autres événements compétitifs dans un modèle classique (Kaplan-Meier, Cox). Avec Fine & Gray, on modèlise tous les risques simultanément. Cela peut arriver notamment lorsque les différents événements compétitifs sont des décès cause-spécifiques : décès pour accident hémorragique, décès pour accident thromobotique, décès pour cancer, etc. Dans ce cas, les événements allant tous dans le même sens (nocifs), l’approche de Kalbfleisch & Prentice est préférable alors qu’avec Fine & Gray, toute exposition qui augmente le risque de décéder d’une cause particulière aurait tendance à donner l’impression qu’il réduit le risque de décès des autres causes.

Synthèse

Quand les événements sont opposés, c’est-à-dire qu’un événement est bénéfique, comme la guérison, alors que l’autre est nocif, comme le décès, alors utilisez l’approche de Fine & Gray pour l’analyse de l’événement bénéfique (guérison). Ce modèle considèrera que les sujets décédés ne guériront jamais. Pour l’analyse du décès, ignorez complètement la guérison, car ce n’est pas un risque compétitif, puisqu’il reste possible d’observer un décès après guérison.

Quand les événements sont tous deux nocifs (p.e. pemphigoïde bulleuse et décès), alors utilisez l’approche de Kalbfleisch & Prentice, c’est-à-dire, analysez le risque de développer une pemphigoïde bulleuse chez les sujets survivants, en censurant sur le décès. Pour cela, vous pouvez utiliser les procédures habituelles des logiciels statistiques (Cox, Kaplan-Meier), en recodant seulement les décès comme des censures de la pemphigoïde bulleuse. Pour l’analyse du décès, ignorez complètement la pemphigoïde bulleuse, car ce n’est pas un risque compétitif, puisqu’il reste possible d’observer un décès après une pemphigoïde bulleuse.

Les flogs : nouveau concept sur un outil ancien

Histoire

John Napier (1550 – 1617) est un mathématicien, physicien et astronome écossais qui inventa le logarithme. La première table de logarithme publiée par Napier était conçue pour effectuer des multiplications par le sinus d’angles et n’utilisait pas la fonction logarithme telle qu’elle est connue de nos jours (référence : Pochon, Luc-Olivier. « A propos d’une table de logarithmes. »).

La fonction logarithme fut extrêmement utile pour simplifier les calculs de multiplication et divisions en de simples additions et soustractions, en s’aidant de tables de logarithme à une époque où les calculatrices n’existaient pas. Ces tables sont restées une manière très efficace et précise de faire des multiplications de nombre décimaux jusqu’au XXième siècle. Les règles à calcul, basées elles aussi sur l’échelle logarithmique permettaient de faire des multiplications et divisions approximatives encore plus rapidement. Ces outils devirent progressivement désuets durant la deuxième moitié du XXième siècle, avec la concurrence des calculatrices mécaniques capables de faire des multiplications, puis dans les années 1970 et 1980, les calculatrices électroniques dont les prix, rapidement décroissants permirent à tous les foyers de s’équiper. Le logarithme et l’exponentielle restent des outils mathématiques indispensables, utilisés au quotidien par le biostatisticien, par exemple pour estimer les paramètres d’une régression logistique ou calculer une variance de Greenwood.

Problématique

Vous avez peut-être déjà atteint les limites des nombres à virgule flottante 64 bits que l’on trouve dans tous les logiciels statistiques. Un nombre inférieur en valeur absolue, à 10^-323 ou supérieur à 10^308, sera représenté comme 0 ou +Inf. Ces limites peuvent paraître, au premier abord, difficiles à atteindre, mais en pratique, on les atteint très vite dès qu’on cumule de nombreuses multiplications ou des exponentiations. Par exemple, dès qu’on calcule la vraisemblance d’un échantillon, on passe rapidement en dessous de 10^-323. Sur un échantillon de taille 1000, avec un pourcentage à 90%, on a une probabilité de 10^-1000 d’observer zéro événement sur les 1000 observations indépendantes, selon la loi binomiale. La probabilité d’observer une valeur au-delà de +100 écarts-types d’une loi normale est d’environ 3×10^-2218. Ce genre de probabilité devra être calculé dans de nombreuses procédures statistiques, telles que les estimateurs du maximum de vraisemblance. C’est une des raisons qui justifie le recours à la log-vraisemblance plutôt que la vraisemblance brute. Grâce à cette transformation logarithmique, les multiplications et les puissances sont respectivement remplacées par des additions et multiplications. On peut alors représenter des vraisemblances de l’ordre de 10^(-10^308).

Conceptualisation de l’outil

Trois points de vue peuvent être posés sur ces transformations logarithmiques :

  1. Il s’agit d’un outil mathématique simplifiant parfois les algorithmes (p.e. dérivation souvent plus simple sur la log-vraisemblance que sur la vraisemblance) ;
  2. Il s’agit d’un outil informatique permettant de dépasser les limitations des nombres à virgule flottante, mais obligeant à réécrire les algorithmes, ce qui est parfois délicat pour les novices ;
  3. Il s’agit d’un détail de représentation interne de l’ordinateur, permettant de gérer des nombres extrêmement grands et petits. Il n’y a pas de mise à jour notable d’un algorithme pour utiliser ces grands nombres.

Les points numéros 1 et 2 correspondent à la vision la plus classique de l’usage du log. Le point numéro 3 est un peu plus original et correspond à un package R que je suis en train de développer.

L’idée, c’est de représenter un nombre quelconque dans l’espace des nombres réels, par deux valeurs:

  1. Une valeur logique (TRUE/FALSE) représentant le signe positif ou négatif du nombre
  2. Une valeur numérique flottante (64 bits) représentant le logarithme du nombre

Une telle représentation sera nommée un flog, mot valise issu de float et log. Ainsi, des nombres aussi petits ou grands que exp(-10^308) et exp(10^308) peuvent être représentés par un flog. Les fonctions arithmétiques de base sont assez facilement réécrites afin de permettre un usage tout à fait normal des flogs, comme s’il s’agissait de valeurs numériques habituelles : additions, soustractions, divisions, multiplications, exponentiations, comparaisons peuvent être réalisées sur les flogs de manière transparente, comme si on manipulait les nombres non transformés. Grâce aux capacités de surcharge des opérateurs de R, un algorithme qui jadis, était programmé pour des nombres à virgule flottante normaux (allant de 10^-323 à 10^308) peut être réutilisé pour les flogs sans modification, gérant ainsi des nombres allant de 10^(-10^308) à 10^(10^308).

Apparement, une recherche Google permet de découvrir que je ne suis pas le premier à avoir eu cette idée, puisqu’un package Haskell existe pour ça (https://hackage.haskell.org/package/logfloat-0.13.4/docs/Data-Number-LogFloat.html).

Outre la simplification donnée par les flogs par rapport à une procédure plus classique de transformation logarithmique, ceux-là ont l’avantage d’être capables de représenter les nombres négatifs comme les nombres positifs.

Rapprochement d’un ancien concept : le float

Après tout, les nombres à virgule flottante (float) ont une représentation interne avec un signe, une mantisse et un exposant (qui se rapproche du log du nombre) de manière transparente. On les manipule comme s’il s’agissait de nombre réels, en ignorant les détails de représentation interne, et notamment, le fait qu’ils sont représentés en base 2. Les flogs poussent le concept un cran au-dessus.

En langage C++, on pourrait écrire un template flog<T>, paramétré par un type numérique, tel que float ou double, qui transformerait ce type numérique en un flog. Comme le flog lui-même est un type numérique, on pourrait créer un flog<flog<double> > qui serait une représentation numérique de très très grand nombres dont le logarithme est lui-même représentable comme un flog. C’est flog<flog<double> > pourraient ainsi représenter des nombres allant jusqu’à 10^(10^(10^308)). On peut appeler ça un flog carré. Bien sûr, on peut aller au-delà, avec les flog<flog<flog<double> > > que l’on appelerait les flogs cubes, et d’une manière plus générale, les flogs itérés à l’ordre k.

Les flogs itérés peuvent paraître superflus, mais mon expérience personnelle m’a montré que l’on pouvait malheureusement dépasser les capacités des flogs simples avec la famille des lois de Weibull. Si on veut représenter correctement la probabilité d’une survie dépassant un certain seuil. Par exemple, un flog est incapable de représenter le taux de survie au temps t=200 d’une distribution de Weibull de paramètre de forme égal à 200 et de paramètre d’échelle égal à 1. Un flog carré est capable de représenter le taux de survie au temps t=10^300 d’une distribution de Weibull de paramètre de forme égal à 10^300. On peut aussi utiliser un cflolog dont je parlerais plus tard.

Limitations générales des floats et flogs

Les limitations des nombres à virgule flottante (floats) sont exacerbés sur les flogs. Tout nombre réel n’est pas représentable avec un float ou un flog. Au maximum 2^64 valeurs distinctes sont représentables. Ainsi, on peut définir la précision d’un nombre float ou flog par la différence qui existe entre les deux nombres successifs représentables sur le float ou flog, sans que ne soit représentable le moindre autre nombre entre les deux. Ainsi, la précision de la valeur 1024 pour un float 64 bits est de 2.3×10^-13, parce que le plus petit nombre strictement supérieur à 1024 représentable par un float 64 bits est égal à (1024+2.3×10^-13). Il s’agit de la précision absolue. La précision relative est obtenue en divisant la précision absolue par le nombre considéré (1024 dans l’exemple). C’est ainsi que le nombre 1024 a une précision relative de 2.2×10^-16.

Contourner les problèmes de précision sur les floats

Les floats ont une précision relative presque indépendante de la valeur; cette précision relative varie entre 1.1×10^-16 et 2.2×10^-16 (sauf pour les dénormaux, mais c’est une autre histoire). On peut donc considérer qu’un float 64 bits a 15 à 16 chiffres numériquement significatifs, peu importe sa valeur. La précision absolue sera très mauvaise pour des valeurs énormes, comme 10^300 qui aura une précision absolue d’environ 10^284, mais sera excellente pour des valeurs petites, comme 10^-280 qui aura une précision absolue d’environ 10^-296. En conséquence, il faut être très délicat avec l’ordre dans lequel on effectue les opérations afin de ne pas perdre en précision inutilement. Imaginons que l’on veuille calculer z=y+x-y/(1+x^2) en sachant que x est petit (p.e. 10^-9) alors que y est grand (p.e. 10^9). Si on calcule cela directement avec des float 64 bits, on obtiendra, à tort, la valeur zéro comme résultat, parce que y+x = y et 1+x^2 = 1 en raison des erreurs d’arrondi. Si on change l’ordre du calcul à z=y-y/(1+x^2)+x, alors on obtiendra, pour résultat, 10^-9, parce qu’on aura supprimé l’erreur d’arrondi correspondant à l’addition de x. Cependant, on n’aura pas supprimé l’erreur d’arrondi conduisant à 1+x^2 = 1 alors qu’en réalité, 1+x^2 = 1+10^-18. Pour corriger cette deuxième erreur d’arrondi, il va falloir réécrire y-y/(1+x^2)+x = y*(1-1/(1+x^2))+x puis se concentrer sur l’estimation de 1-1/(1+x^2). Considérant que x^2 est très proche de zéro, on peut approximer la fonction inverse à son développement limité à l’ordre 1, et donc, approximer 1/(1+x^2) à 1-x^2; d’où une approximation de 1-1/(1+x^2) à x^2. On réécrit alors z ≃ y×x^2+x qui se calcule alors de manière précise à 2×10^-9. Paradoxalement, en utilisant un algorithme approximatif (développement limité), on obtient un résultat plus précis qu’avec un algorithme « exact » à ceci-près que ce dernier ignore le problème d’imprécision des nombres à virgule flottante. En pratique, si on veut implémenter cela de manière générique et propre, on utilisera un développement limité à un plus grand ordre (p.e. 1-1/(1+x^2) ≃ x^2 – x^4 +x^6 – x^8) et on définira un seuil sur la valeur de x (p.e. 10^-2) en-dessous duquel on utilise le développement limité, et au-dessus duquel, on utilise le calcul direct 1-1/(1+x^2). Ce seuil doit être choisi de telle sorte que la précision des deux méthodes soit satisfaisante à ce seuil. Ici, l’erreur absolue du développement limité est environ égal à x^10, c’est-à-dire, 10^-20 au niveau du seuil, soit une erreur relative d’environ 10^-16. C’est donc une précision satisfaisante. Pour le rapport 1-1/(1+x^2) calculé directement, l’erreur absolue de 1+x^2 est d’environ 10^-16, ce qui conduit à une erreur absolue du même ordre de grandeur pour son inverse ainsi que pour 1-1/(1+x^2) et donc, une erreur relative d’environ 10^-12. La précision est moins bonne, mais reste satisfaisante. Les développements limités s’avèrent très précieux dans la manipulation des nombres à virgule flottante. Pour cela, il ne faut pas hésiter à s’aider du logiciel libre Xcas qui calcule assez bien les développements limités de manière formelle. Pour certaines fonctions fréquemment utilisées, il existe des procédures dans la plupart des logiciels et langages de programmation : expm1 pour exp(x)-1, log1p pour log(1+x). De manière un peu plus rare, on trouve log1pmx(x)=log(1+x)-x et lgamma1p(x)=log(gamma(1+x)). Ces fonctions sont faciles à implémenter avec des développements limités si elles ne sont pas nativement disponibles dans le logiciel.

Limitations supplémentaires des flogs

La précision absolue des floats est décroissante avec leur taille alors que la précision relative est à peu près constante, aux alentours de 2×10^-16 pour les floats 64 bits. Avec les flogs, la précision relative devient décroissante avec la taille des nombres. Ainsi, la table suivante donne des points de repères :

ValeurPrécision relative d’un float 64 bitsPrécision relative d’un flog 64 bits
10^-(10^100)Non représentable*10^(3.9×10^84)
10^-(10^16)Non représentable*54
10^-(10^15)Non représentable*0.65
10^-(10^10)Non représentable*3.8×10^-6
10^-1000Non représentable*4.5×10^-13
10^-100~ 2.2×10^-162.8×10^-14
10^-15~ 2.2×10^-167.1×10^-15
1~ 2.2×10^-164.9×10^-324
10^15~ 2.2×10^-167.1×10^-15
10^100~ 2.2×10^-162.8×10^-14
10^1000Non représentable**4.5×10^-13
10^(10^10)Non représentable** 3.8×10^-6
10^(10^15)Non représentable** 0.65
10^(10^16) Non représentable** 54
10^(10^100)Non représentable** 1.94×10^84
Précisions relatives des floats et flogs 64 bits

* Non représentable ou représenté comme strictement égal à 0, et donc, infiniement imprécis, selon le point de vue.

** Non représentable ou représenté comme égal à +Inf, et donc, infiniement imprécis, selon le point de vue.

La précision relative des flogs est optimale autour de 1, bien meilleure que celle des floats. La précision relative des flogs reste excellente jusqu’à 10^1000 et est tolérable jusqu’à 10^(10^10). À partir de 10^(10^16), la précision relative descent à 54, ce qui veut dire que le rapport entre deux flogs successifs est égal à 55. En conséquence, en considérant un mode de calcul qui arrondit au nombre représentable le plus proche, si on additionne un nombre avec lui-même, il est inchangé (x+x = x), et si on multiplie un nombre de l’ordre de 10^(10^16), par un facteur inférieur à 27.5, il est inchangé (p.e. x×20 = x). En première impression, l’addition n’a plus aucun sens sur des valeurs aussi élevées, puisque l’addition de deux nombres de cet ordre de grandeur, sera toujours égal au plus grand des deux. Pourtant, il est m’est arrivé d’additionner et soustraire des flogs bien plus grands, sans que cela pose problème. Lorsqu’on veut savoir si 10^(10^270)+ 10^(10^300) est supérieur ou inférieur à 10^(10^271)+10^(10^301), on comprend que les « petits » nombres 10^(10^270) et 10^(10^271) sont totalement négligeables dans la comparaison et peuvent être ignorés.

Pour 10^(10^100), on voit que la précision relative est très mauvaise, puisque deux nombres représentables successifs ont un rapport 1.94×10^84 l’un avec l’autre. Si on multiplie un nombre de cette taille par des milliards de milliards, il reste inchangé (x×10^18 = x). Pourtant, certaines opérations restent signifiantes à cette échelle : un nombre élevé au carré sera toujours estimé avec une bonne précision. En fait, pour passer d’un nombre représentable au plus petit nombre représentable strictement supérieur, il faut élever le premier à la puissance (1+2.2×10^-16) environ. L’opérateur d’exponentiation (^) garde alors pleinement son sens.

Exemple d’implémentation de l’addition de flogs

Comment peut-on calculer a+b à partir de la=log(a) et de lb=log(b) et le représenter sur une échelle logarithmique ? Si on calcule directement log(exp(la)+exp(lb)), les représentations intermédiaires de exp(la) et exp(lb) dépassent les capacités des floats. Il existe pourtant une stratégie simple et efficace de calcul, basé sur le fait que pour toute constante c > 0, exp(la)+exp(lb) = (exp(la)/c + exp(lb)/c)×c et donc log(exp(la)+exp(lb)) = log(exp(la – log(c)) + exp(lb – log(c))) + log(c). Il faut ensuite choisir une constante lc = log(c) qui garantisse que l’on ne dépasse pas les capacités des nombre à virgule flottante. Pour cela, la constante lc=max(la, lb) est idéale. Ainsi, la-lc et lb-lc seront tous deux inférieurs ou égaux à 0 et on ne dépassera pas les capacités des nombre à virgule flottante. Si on additionne des nombres qui n’ont pas le même ordre de grandeur, tels que e^10000 et e^10100, le plus petit des deux sera négligeable face au plus grand, et le résultat final sera égal au plus grand des deux, en raison des arrondis flottant.

Place des flogs et floats

Les floats sont tout à fait adaptés à des algorithmes ayant un recours massif aux additions et comprenant bon nombre de multiplications, mais on arrive à leurs limites dès que l’on commence à utiliser des puissances ou que le recours aux multiplications est trop massif. C’est là que les flogs ont leur place.

Bon usage des flogs

Les limites des flogs sont semblables à celles de floats, mais il est encore plus nécessaire d’en avoir conscience car ils sont conçus pour manipuler des valeurs extrêmes. C’est ainsi que contrairement à ce que j’ai décrit en introduction, il faut avoir conscience de la représentation des nombres pour écrire des algorithmes corrects. Par exemple, si x est potentiellement immensément grand, positif ou négatif, il faudra absolument éviter d’écrire des expressions telles que log(1+exp(x)), puisque cette expression dépassera les capacités des nombres pour un grand x positif et perdra toute précision pour un grand x négatif. Pour x positif et supérieur à 37 environ, on pourra approximer log(1+exp(x)) à x, sans perte de précision. Pour x négatif et inférieur à 37 environ, on pourra approximer log(1+exp(x)) à exp(x). Pour x entre les deux, on pourra s’aider de la fonction log1p pour écrire log1p(exp(x)).

Pour aller plus loin

La transformation logarithmique n’est pas la seule qui soit possible pour créer un nouvel espace de nombres à partir des floats. On peut notamment, pour des valeurs comprises entre 0 et 1 utiliser la transformation cloglog qui à x associe log(-log(x)) pour construire ainsi les cflologs. Seuls les nombres compris entre 0 et 1 sont représentables ainsi.

Là où les floats sont incapables de représenter des nombres compris entre 1-1.1×10^-16 et 1, les flogs sont capables de représenter 1-5×10^-325 et les cflologs sont capables de représenter 1-10^(-1.8×10^308). Là où le plus petit nombre strictement positif représentable par les floats est 5×10^-325, les flogs peuvent représenter 10^(-1.8×10^308) et les cflologs peuvent représenter 10^(-10^(1.8×10^308)). On peut en effet dépasser les capacités des flogs avec la distribution de Weibull, parce que celle-ci dépend un exponentiation dans une exponentielle, la fonction cumulative de distribution de Weibull étant F(x; forme, échelle) = 1-exp(-(x/échelle)^forme).

Le problème de cette transformation, c’est que l’espace de représentations numériques ainsi créé est limité à l’intervalle [0, 1], interdisant alors bon nombre d’opérations. Si on veut rendre ces représentations polyvalentes, il faut enrichir la représentation. Par exemple, on pourrait proposer deux floats x et y, représentant la valeur exp(-exp(x))×y.

Conclusion

Les flogs sont un nouveau point de vue, sur un outil ancien. Utilisez-les à bon escient !

Que faire d’un paramètre de nuisance ?

Le problème d’estimation d’un simple pourcentage est relativement simple. En effet, on dispose généralement de deux données : un numérateur et un dénominateur du pourcentage. La distribution du numérateur suit une loi binomiale qui est parfaitement déterminée pour toute proportion théorique définie. Les fluctuations d’échantillonnages sont ainsi facilement identifiables et, par inversion de tests binomiaux, on peut aisément construire un intervalle de confiance d’un pourcentage. Il y a des subtilités liées à l’aspect discret de la loi binomiale qui rend impossible le calcul d’un intervalle de confiance exact, mais globalement, ça reste encore simple.

Le problème d’estimation d’une différence de deux pourcentages est nettement plus complexe car les fluctuations d’échantillonnages de la différence des deux pourcentages ne peut pas être déterminé par cette différence seule. Par exemple une différence de 10% entre deux pourcentages peut exister pour deux proportions théoriques π1=0.10 et π2=0.20 aussi bien que pour π1=0.50 et π2=0.60. Dans le premier cas, l’estimateur empirique de la différence aura une variance beaucoup plus faible que dans le second cas, car les deux proportions théoriques ont chacune une variance moindre. Lorsque la distribution du paramètre que l’on souhaite estimer n’est pas seulement dépendante du paramètre lui-même mais dépend d’autres paramètres, inconnus, ces autres paramètres inconnus sont nommés paramètres de nuisance. Comment peut-on alors évaluer l’incertitude sur l’estimation d’un paramètre en présence d’un ou plusieurs paramètres de nuisance ?

Deux stratégies classiques existent, mais j’utilise parfois une troisième :

  1. Le conditionnement sur le paramètre de nuisance observé
  2. L’usage de la pire valeur théoriquement possible du paramètre de nuisance
  3. Une approche mixte bayésienne-fréquentiste prenant en compte la distribution bayésienne a posteriori du paramètre de nuisance, compte-tenu de l’échantillon, pour un prior non informatif.

Exemple

L’exemple de l’estimation de la différence de deux pourcentages sera utilisée pour la suite de ce billet de blog. On note n1 et n2 les tailles de deux échantillons (dénominateurs des pourcentages). On note x1 et x2 les numérateurs observés des pourcentages. On note π1 et π2 les vraies proportions dans la population alors que p1=x1/n1 et p2=x2/n2 sont les proportions observées sur les échantillons. Les échantillons sont supposés être constitués par tirage au sort d’observations indépendantes et identiquement distribuées dans chacun des deux échantillons.

On souhaite estimer π1 – π2 et le paramètre de nuisance est (π1 + π2)/2, la moyenne des deux vraies proportions. Enfin, on note X1, X2, P1 et P2 les variables aléatoires dont les réalisations sont x1, x2, p1 et p2.

Conditionnement

Cette approche consiste à supposer que le paramètre de nuisance réel est égal au paramètre de nuisance estimé sur l’échantillon. Dans l’exemple, on va donc supposer que (π1 + π2)/2 est strictement égal à (p1 +p2)/2 sans incertitude.

Ceci est clairement une approximation, dont les conséquences vont être un biais statistique plus ou moins important.

Toujours sur le même exemple, la variance de la différence des deux pourcentages sera maximale dans le cas où π1 et π2 sont tous deux égaux à 0.50. Ainsi, d’une manière générale, la variance de P1-P2, sera plus grande pour (π1 + π2)/2=0.50 que pour les autres valeurs. Or, même si (π1 + π2)/2=0,50, la valeur de (p1 + p2)/2 sera différente et correspondra à une variance plus faible de P1-P2 . C’est ainsi qu’on sous-estimera systématiquement la variance de P1-P2; cela biaisera les intervalles de confiance et petits p.

Dans la plupart des cas, les sous-estimations ou surestimations du paramètre de nuisance ont tendance à compenser les biais, en moyenne, et le conditionnement sur le paramètre de nuisance n’est pas forcément associé à un gros biais. Les choses sont généralement plutôt biaisées lorsqu’il existe une covariance importante entre le paramètre que l’on souhaite estimer et le paramètre de nuisance.

La méthode du conditionnement est extrêmement fréquemment utilisée, valable sur des échantillons moyens ou grands, mais souvent biaisée sur de petits ou très petits échantillons.

On notera que certaines méthodes se vantant d’être « exactes » comme la méthode de régression logistique exacte d’Hirji (https://doi.org/10.2307/2289388) s’avèrent en réalité être conditionnelles aux paramètres de nuisance; les autres coefficients du modèle pour la méthode d’Hirji.

Pire valeur théoriquement possible

Un certain nombre de méthodes se vantant d’être « exactes » résolvent le problème du paramètre de nuisance en se basant sur la valeur du paramètre pour laquelle l’incertitude du paramètre estimé est la plus grande. Avec cette stratégie du « pire des cas », la méthode est forcément conservative, c’est-à-dire, que les niveaux de confiance réels des intervalles de confiance sont supérieurs ou égaux au niveau nominal et les risques alpha réels sont inférieurs au risque nominal.

L’archétype de cette stratégie est illustrée par la méthode exacte de Santner et Snell (https://doi.org/10.1080/01621459.1980.10477482) pour l’estimation de la différence de deux pourcentages ou le rapport de deux pourcentages. Cette méthode est un bon exemple des effets pervers de cette stratégie. Les intervalles de confiance de Santner et Snell sont extraordinairement conservateurs. Quand je dis extraordinairement, c’est, par exemple, que le rapport de la proportion 20/40 avec la proportion 20/40, fournit l’intervalle de confiance à 95% [0 à +infini]. Je n’ai pas vérifié avec un logiciel (en raison du temps de calcul trop long), mais je crois que même avec 50000/100000 dans chacun des deux échantillons, l’intervalle de confiance serait [0 à +infini].

Pourquoi la méthode de Santner & Snell en arrive-t-elle a une telle aberration ? C’est parce qu’elle considère que même si le paramètre de nuisance observé est de 50% pour ce cas là, le paramètre de nuisance théorique le pire est à peu près à (π1 + π2)/2 =1/40. En effet, pour une telle valeur du paramètre de nuisance, les fluctuations d’échantillonnage sont absolument énormes puisque le nombre d’événements attendu dans chacun des groupes est égal à 1, avec une haute probabilité de division par zéro, qui conduit à un risque relatif infini. Le problème de l’approche de Santner & Snell, c’est qu’elle va s’intéresser à un paramètre de nuisance totalement extrême mais pas du tout plausible compte-tenu de l’échantillon. Cela conduit à un conservatisme extrême qui rend la méthode totalement inutile et non pertinente.

La plupart du temps, on ne tombe pas dans ces extrémités, mais je voudrais néanmoins argumenter en défaveur des approches conservatrices. En effet, un fréquentiste conservateur aura tendance à considérer que, si le niveau de confiance réel d’un intervalle de confiance peut, selon un paramètre inconnu, valoir 98% ou 92%, alors on peut considérer qu’il s’agit d’un intervalle de confiance à 92%. De mon point de vue, si on peut raisonnablement parier (point de vue bayésien) qu’il y a une chance sur deux que le niveau de confiance soit de 98% et une chance sur deux que le niveau de confiance soit 92%, alors le niveau de confiance est de 95%. Pour caricaturer encore le propos, imaginons qu’on tire à pile ou face pour déterminer si on va calculer un intervalle de confiance à 92% ou à 98%. La procédure complète, du tirage au sort, suivi du calcul, correspond, de mon point de vue, à la construction d’un intervalle de confiance à 95%.

Cela permet d’introduire la troisième approche, qui vise justement à prendre en compte la distribution bayésienne du paramètre de nuisance.

Approche mixte bayésienne-fréquentiste

Cette approche consiste à considérer l’ensemble des valeurs plausibles du paramètre de nuisance, associées à leurs probabilités. Même si cela combine de manière non rigoureuse la théorie fréquentiste et bayésienne, les méthodes produites sont tout à fait analysables et interprétables de manière fréquentistes pures; on peut dire que l’aspect bayésien n’est qu’un détail de construction. La procédure est la suivante:

  1. Estimer la distribution bayésienne du paramètre de nuisance, compte-tenu de l’information de l’échantillon, pour un prior non informatif; cela peut se faire de manière purement fréquentiste, par l’usage d’une distribution de confiance.
  2. Approximer cette distribution du paramètre de nuisance éventuellement continue, à une distribution discrète de valeurs très rapprochées et nombreuses : par exemple, par une distribution discrète à 1024 valeurs.
  3. Pour chacune des 1024 valeurs possibles du paramètre de nuisance, estimer le paramètre d’intérêt, conditionnellement à cette valeur du paramètre de nuisance, soit par une distribution de confiance, soit, si on est paresseux, par les bornes d’un intervalle de confiance à 95%
  4. Calculer la moyenne des 1024 estimations du paramètre d’intérêt pondérée par la probabilité de présence de la valeur du paramètre de nuisance

Cette approche est bien plus fine que la méthode du conditionnement car elle prend en compte l’incertitude sur l’estimation du paramètre de nuisance, avec une probabilité maximale à proximité de l’estimation ponctuelle. Le conditionnement considère que le paramètre de nuisance réel est égal à son estimation ponctuelle sans la moindre incertitude. L’usage de la pire valeur, quant à elle, se base sur des valeurs non plausibles du paramètre de nuisance, et est excessivement conservatrice.

Cette approche a quand même quelques défauts : elle est plus élaborée que le conditionnement, est beaucoup plus lourde en temps de calcul sur ordinateur, puisqu’on va typiquement multiplier les par 1024 par rapport à une approche purement conditionnelle. Elle fournit néanmoins des « diagnostics » intéressants parce qu’elle permet de voir à quel point, sur l’espace des valeurs plausibles du paramètre de nuisance, la valeur de ce paramètre de nuisance influence les résultats. Cela aide à voir à quel point le paramètre de nuisance perturbe les résultats et à quel point on peut craindre que les résultats sont biaisés. Si tout se passe bien, on verra que l’estimation du paramètre d’intérêt est faiblement dépendant du paramètre de nuisance, sur la zone plausible de ce paramètre. Néanmoins, ce type d’approche reste intéressant lorsque la dépendance est assez forte, puisqu’on peut espérer qu’elle soit moins biaisée que le conditionnement classique.

Il est aussi souhaitable, si on développe ce type d’estimateur, de le valider à partir de simulations, idéalement les plus proches possibles du cas de figure auquel on fait face.

Conclusion

WordPress m’a fait perdre 30 minutes de travail en bousillant une sauvegarde de ce billet de blog; c’est vraiment ennuyeux de devoir retaper de mémoire ce qu’on avait écrit la première fois.

Au total, je considère que l’approche mixte bayésienne-fréquentiste est théoriquement intéressante et très générale. Elle me paraît particulièrement pertinente sur de petits échantillons lorsqu’on craint qu’un paramètre de nuisance biaise fortement les statistiques. Elle peut parfois être un peu délicate à mettre en oeuvre, mais c’est aussi un exercice intéressant.

Enfin, cerise sur le gâteau, avec ce type d’approche, vous allez vous faire de nombreux ennemis dans les deux principales religions statistiques : fréquentiste et bayésienne. Avec un peu de chance, au 22ème siècle, les églises seront unifiées et plus personne ne se tapera sur la figure. Peut-être aurez-vous un peu contribué, en votre temps, à la grande unification; même si le martyre reste l’événement le plus probable.

Le code source de R

Je commence à connaître un peu le code source de R, en raison d’améliorations que j’y ai apportées : possibilité d’interrompre les gros calculs, accélération de la génération de nombre aléatoires, suppression de la limite à 128 connections, accélération de la fonction mean() sur les ALTREPs, refactoring de la fonction sum().

Voilà mon opinion : R contient des algorithmes conçus pour une précision numérique maximale, souvent au détriment des performances, ce qui se comprend étant donné les objectifs du projet. Par exemple, dès que possible, des nombres à virgule flottante 80 bits sont utilisés dans les calculs intermédiaires, quand bien même le résultat final est seulement stocké sur 64 bits. Cependant, la structure du code source est vraiment mauvaise. Cela rend le code moins maintenable et plus obscur. Voilà, les éléments que j’ai constatés en lisant le code source des quelques fichiers dont j’ai eu besoin pour mes patches:

Un peu de Fortran 77

Pour commencer par un des aspects les plus acceptables du code source, celui-ci contient du code Fotran 77 avec le style de l’époque, notamment le respect de la limite à 6 caractères pour la longueur des identifiants, y compris des variables locales. Ce n’est pas très agréable à lire, mais heureusement, cela ne correspond qu’à une minorité du code source, et le code est basé sur des interfaces déjà connues (telles que BLAS/LAPACK, LINPACK), même si cela a engendré des problèmes d’interopérabilité Fortran/C.

Il contient aussi un peu de Fortran 90 qui le rend impossible à compiler, pour l’instant, sur Apple M1; d’ici quelques mois il est probable que GCC soit porté sur cette nouvelle plateforme.

Des gros switch() partout

Dans le fichier src/main/RNG.c contenant les 8 algorithmes de génération de nombres pseudo-aléatoires (dont l’algorithme « user-supplied »), les algorithmes sont implémentés de manière éparse dans diverses fonctions.

C’est-à-dire qu’un algorithme devrait implémenter:

  1. Une fonction d’initialisation de l’état du générateur à partir d’une graine pseudo-aléatoire 32 bits (RNG_Init)
  2. Une fonction de remise en cohérence de l’état du générateur au cas où l’utilisateur aurait modifié l’état interne par un accès direct à l’objet .Random.seed (FixupSeeds).
  3. Une fonction de génération de nombre pseudo-aléatoire proprement dite (unif_rand)
  4. Des informations telles que le nombre de bits de précision effective du nombre pseudo-aléatoire

Plutôt que de définir une interface, sous forme d’une structure avec pointeurs vers fonction, l’ensemble des fonctionnalités est implémentée sous forme de gros switch() pour chacune des méthodes.

double unif_rand(void)
{
    double value;

    switch(RNG_kind) {

    case WICHMANN_HILL:
        I1 = I1 * 171 % 30269;
        I2 = I2 * 172 % 30307;
        I3 = I3 * 170 % 30323;
        value = I1 / 30269.0 + I2 / 30307.0 + I3 / 30323.0;
        return fixup(value - (int) value);/* in [0,1) */

    case MARSAGLIA_MULTICARRY:/* 0177777(octal) == 65535(decimal)*/
        I1= 36969*(I1 & 0177777) + (I1>>16);
        I2= 18000*(I2 & 0177777) + (I2>>16);
        return fixup(((I1 << 16)^(I2 & 0177777)) * i2_32m1); /* in [0,1) */

/* ... */

Cela rend l’ajout d’algorithmes très délicat, parce qu’il faut rajouter une entrée dans chacun des switch() de src/main/RNG.c.

Par ailleurs, la liste des générateurs pseudo-aléatoires est identifié par une énumération:

typedef enum {
    WICHMANN_HILL,
    MARSAGLIA_MULTICARRY,
    SUPER_DUPER,
    MERSENNE_TWISTER,
    KNUTH_TAOCP,
    USER_UNIF,
    KNUTH_TAOCP2,
    LECUYER_CMRG
} RNGtype;

Mais dans la fonction PutRNGstate, il est vérifié que le numéro du générateur est « valide »; cette vérification est simplement la comparaison de RNG_king à la valeur LECUYER_CMRG. En bref, le fait que l’algorithme LECUYER_CMRG est le dernier de la liste est hard-codé dans le code source, obligeant à modifier cette partie du code source lors de l’ajout d’un nouveau générateur pseudo-aléatoire.

Des constantes magiques

Les fonctions R de base (primitives) mean, sum, prod, min et max sont implémentées dans la même fonction C : do_summary dans src/main/summary.c. Les cinq fonctions (mean, sum, prod, min et max) sont distinguées par un numéro d’opération (0=sum, 1=mean, 2=min, 3=max, 4=prod) et des gros switch() disséminés dans la fonction do_summary adaptant le comportement à l’opération. Plutôt que de définir une énumération permettant d’associer un numéro à un nom (p.e. MIN_IOP = 2), le code contient des constantes magiques, 0, 1, 2, 3 et 4, conduisant à ce genre de bout de code :

    if(empty && (iop == 2 || iop == 3)) {
	if(ans_type == STRSXP) {
	    warningcall(call, _("no non-missing arguments, returning NA"));
	} else {
	    if(iop == 2)
		warningcall(call, _("no non-missing arguments to min; returning Inf"));
	    else
		warningcall(call, _("no non-missing arguments to max; returning -Inf"));
	    ans_type = REALSXP;
	}
    }

Si on ne sait pas que 2=min et 3=max, on peut difficilement comprendre ce code.

Dans la même fonction do_summary, il y a une variable « empty », partagée entre toutes les opérations, valant TRUE tant qu’on n’a rencontré que des vecteurs vides (longueur zéro) ou ne présentant que des données manquantes (NA) éliminées par un na.rm=TRUE. Cela permet aux fonctions min et max d’émettre un message d’avertissement dans ce genre de cas:

> min(c(NA,NA),NULL,na.rm=TRUE)
[1] Inf
Warning message:
In min(c(NA, NA), NULL, na.rm = TRUE) :
  no non-missing arguments to min; returning Inf

Lors du traitement de chaque paramètre de la fonction, une variable « updated » est mise à jour afin de dire si des valeurs non manquantes (sauf si na.rm=FALSE) ont été intégrées dans l’opération. Eh bien, dans certaines situations la variable updated, qui est booléenne à la base, peut devenir un trooléen pour décrire une situation de dépassement de capacité (overflow) des nombres entiers 64 bits, concernant la somme de nombres entiers. La troisième valeur, qui n’est ni TRUE, ni FALSE, vaut 42, sans la moindre explication quant à ce choix, même si on devine qu’il y a un rapport direct à la grande question sur la vie, l’univers et le reste.

Du code dupliqué en fonctionnalité

Résumé

R contient deux implémentations de fonctions associées aux distributions aléatoires (rnorm, rbinom, …), différant de manière subtile et accidentelle en comportement. Une des deux implémentations n’est jamais utilisée, sauf dans le cas où l’utilisateur outrepasserait volontairement les fonctions proposées par R (p.e. stats::rnorm), afin d’accéder directement aux fonctionnalités internes (.Internal(rnorm(a,b,c)).

Détail

L’implémentation des fonctions de hasard est séparée en deux parties :

  1. D’un côté, les algorithmes pseudo-aléatoires élémentaires, implémentés dans src/main/RNG.c : ceux-ci permettent de générer des nombres selon une loi uniforme entre 0 et 1 (à quelques détails près) ;
  2. D’un autre côté, des algorithmes permettant de générer des distributions particulières (p.e. binomiale, Student, exponentielle, etc.), conçues pour être directement accessibles depuis du code R (rbinom, rt, rexp, etc.).

En apparence, la deuxième partie est implémentée dans le fichier src/main/random.c, sauf que la grande majorité du code de ce fichier ne semble servir à rien ! Il ne semble jamais être appelé, jamais être utilisé ! La seule partie utilisée concerne la fonction R sample(), mais tout ce qui concerne rchisq, rexp, rgeom, rpois, rt, rsignrank, rbeta, rbinom, rcauchy, rf, rgamma, rlnorm, rlogis, rnbinom, rnorm, runif, rweibull, rwilcox, rhyper, n’est jamais appelé par du code R (nous verrons par la suite quand ça peut encore être utilisé). J’ai pu supprimer 257 lignes de code de src/main/random.c, sans perte évidente de fonctionnalité de R. Comment cela est-il possible ? Eh bien, ces fonctionnalités ont été réimplémentées dans src/library/stats/src/random.c, compilé dans une bibliothèque dynamique stats.so ou stats.dll qui sera ensuite chargée en même temps que le packages ‘stats’ qui fait partie des package pré-chargés de R. Cette réimplémentation est quasiment identique, en fonctionnalité, à l’implémentation de src/main/random.c, avec néanmoins de très subtiles différences, telles que la fonction rbinom() qui renverra des nombres entiers ou flottants, selon qu’ils dépassent ou pas INT_MAX (2147483647) dans src/library/stats/src/random.c, alors qu’elle renverra toujours des flottants dans src/main/random.c. Même si les deux implémentations sont presque identiques en fonctionnalité, elles diffèrent nettement en code source, avec néanmoins des points de convergence assez nets.

On peut en comprendre l’histoire, en remontant l’historique des commits. En avril 2012, il n’y avait pas de doublon de fonctionnalité, puisque src/library/stats/src/random.c ne contenait que le code de rmultinom et r2dtable qui ne sont pas implémentés dans src/main/random.c. En mai 2012, le code de src/main/random.c a été dupliqué dans src/library/stats/src/random.c. Déjà, à l’époque, la duplication du code n’était pas fidèle : c’était une réimplémentation compatible mais de style différent. Il y avait déjà des divergences de fonctionnalité apparente; notamment rbinom() renvoyait toujours des nombres entiers avec l’implémentation de src/library/stats/src/random.c, avec possible overflow, alors que le code de src/main/random.c renvoyait des nombres à virgule flottante. Par la suite les deux implémentations ont divergé de plus en plus, y compris avec des patches qui modifiaient le code de src/main/random.c sans savoir qu’il n’était pas vraiment utilisé. Le grand responsable semble donc être le patch de mai 2012. En même temps que de dupliquer tous les algorithmes de génération de nombre aléatoires (rbinom, rchisq, …), ce patch a aussi dupliqué tous les algorithmes de densité de probabilité, de fonction de répartition et de fonctions des quantiles (dnorm, pnorm, qnorm, dbinom, dbinom, …) en créant le fichier src/library/stats/src/distn.c qui faisait encore 481 lignes de fonctionnalité dupliquée. Ce patch semble correspondre au CHANGELOG de R 3.0.0 qui précise (citation exacte):

CODE MIGRATION

  • The C code underlying base graphics has been migrated to the graphics package (and hence no longer uses .Internal() calls).
  • Most of the .Internal() calls used in the stats package have been migrated to C code in that package. This means that a number of .Internal() calls which have been used by packages no longer exist, including .Internal(cor) .Internal(cov), .Internal(optimhess) and .Internal(update.formula).
  • Some .External() calls to the base package (really to the R executable or shared library) have been moved to more appropriate packages. Packages should not have been using such calls, but some did (mainly those used by integrate()).

Les fonctions .Internal sont celles que l’on retrouve dans le dossier src/main, telles que les fonctions de src/main/random.c. Cela suggère donc une migration de fonctionnalité dont l’intérêt n’est pas complètement évident, mais cela n’explique pas la duplication de fonctionnalité, qui a mené à de dangereuses divergences de code par la suite. Pourquoi cette duplication ? On a un élément de réponse dans le message venant avec le commit de mai 2012:

port C code for [dpqr]xxx functions to stats
leave .Internals for now, as they are used directly as an optimization by byte-compiler

Ainsi, le code inutilisé en apparence pourrait encore être utilisé par le byte-compiler; sauf que c’est faux, bien heureusement. Le byte-compiler est une fonction introduite dans R 2.13.0, avec la précompilation de tous les packages par défaut à partir de R 2.14.0. C’est R 2.15.0 et R 3.0.0 que le patch de mai 2012 a été introduit. Cette compilation convertit le code R en code de machine virtuelle basée sur une pile (stack based virtual machine), avec une amélioration très légère en performances, notamment en raison de certaines optimisations. Ces optimisations consistent en le remplacement de certaines fonctions de bases (p.e. addition, boucle for, parenthèse ouvrante, accolade ouvrante) par une référence directe à la fonction interne (builtin ou special), outrepassant les accès aux environnements de travail, qui sont très chronophages. Dans l’idée du commentaire du commit de mai 2012, le byte-code permettrait l’accès direct aux fonctions telles que rnorm() sans passer par la version réécrite pour le package stats.

Mais cela est presque totalement faux, pour plusieurs raisons :

  1. Le code responsable de la compilation dans src/library/compiler/R/cmp.R n’outrepassera le code de la fonction stats::rnorm() et autres fonctions aléatoires seulement dans le cas où le corps de celle-ci serait limité à un appel direct à .Internal, grace à la vérification faire dans is.simpleInternal. Ainsi, le byte-code compilé par R 3.0.0 ne fera pas d’optimisation et appellera toujours la version implémentée dans src/library/stats/src/random.c.
  2. R 3.0.0 refusera de charger un package compilé pour R 2.15, évitant ainsi tout risque d’utiliser du byte-code optimisé pour la version précédente et faisant un accès direct à la version interne implémentée dans src/main/random.c.

Il reste théoriquement possible d’exécuter l’ancienne version avec la procédure suivante

  1. Appeler compiler::setCompilerOptions(optimize=3) sous R 2.15
  2. Compiler du byte-code toujours sous R 2.15 et l’exporter dans un fichier rds avec saveRDS()
  3. Charger ce byte-code sous R 3.0.0 avec readerRDS()
  4. Évaluer ce code avec eval()

Compiler du code sous une version de R et l’exécuter sous la suivante est à coup sûr une mauvaise idée de toute façon. Évidemment, il existe aussi la possibilité qu’un programme accède directement et explicitement à la fonction interne avec un appel à .Internal. Il n’est pas étonnant que ce genre de programme se casse dans une nouvelle version de R. En pratique, le byte-code compilé avec R 2.15 s’exécutera correctement sur une version récente de R (p.e. 4.1.0) parce que le byte-code ayant évolué entre temps, R ignorera le byte-code de l’expression compilée et se contentera d’utiliser l’expression brute, non compilée, évitant alors d’exécuter la version de src/main/random.c.

Dans tout les cas, la solution retenue est la pire qui soit. Dans les très rares cas où la fonction interne de src/main/random.c serait appelée, par l’exécution de code byte-compilé pour la version précédente de R, ou à cause d’un appel direct par .Internal, le comportement va diverger de celui de la fonction standard (maintenant src/library/stats/src/random.c) en exécutant un code C de plus en plus mal testé au gré des versions successives de R.

Les solutions pertinentes auraient été:

  1. Supprimer le code de src/main/random.c et les fonctions internes correspondantes afin de ne conserver que la nouvelle implémentation avec src/library/stats/src/random.c
  2. Remplacer le code de src/main/random.c par des redirections vers les implémentations de src/library/stats/src/random.c conservant ainsi la compatibilité avec les cas d’usage les plus bizarres
  3. Ne pas migrer le code dans src/library/stats/src/random.c

La même réflexion s’applique aux fonctions dupliquées dans src/library/stats/src/distn.c et src/main/arithmetic.c.

On peut aussi s’étonner de l’incohérence avec la migration d’autres parties du code, telles que les fonctions cor et cov dont l’implémentation interne dans src/main avait été totalement supprimée lors de la migration. La raison semble être due à la tentative d’optimisation du byte-code qui n’existait pas pour ces dernières fonctions.

J’ai mentionné ce problème sur la mailing list r-devel, mais ça ne semble pas intéresser les développeurs R, le message étant resté sans réponse.

Conclusion

Le code de R aurait besoin d’un refactoring important, mais ce n’est pas le style de la maison. Les recommandations de développement de R sont claires :

DO NOT fix exotic bugs that haven't bugged anyone
DO NOT fix things that are not broken
DO NOT restructure code

Mise à jour

J’ai rapporté le bug du code dupliqué sur le bugzilla de R et ai rapidement reçu la réponse des développeurs de R:

WONTFIX « I prefer to leave this as is for now.« 

Ce qui signifie que l’équipe de développement a conscience de ce problème mais préfère conserver le code dupliqué, tel qu’il est; probablement par peur de créer un problème de compatibilité, sans bénéfice immédiat évident. Cela est en cohérence avec la politique :

DO NOT fix things that are not broken
DO NOT restructure code

Adoption jeux d’instructions

Un petit billet sur un sujet informatique pour changer des sujets plus statistiques. J’ai fait un récapitulatif de la gestion des jeux d’instructions complémentaires (SSE2, SSSE3, SSE4.2, AVX, AVX2, AVX-512) par les processeurs AMD et Intel. Cette liste devrait permettre aux programmeurs de savoir ce qu’on peut raisonnablement considérer comme universellement présent et ce qui est seulement disponible sur certains processeurs uniquement. Par exemple, il est possible que le jeu vidéo Star Citizen requière AVX par méconnaissance du défaut de support par les processeurs Pentium Gold récents; ce n’est pas le haut de gamme, mais couplé à un GPU correct, cela aurait dû faire tourner tous les jeux vidéos.

Tableaux bruts

Jeu d’instructionsFamilleDate d’introductionDernier proc à ne pas gérerDébut du support complet sur toute la gamme
SSE2Atom2008 (first Atom)Aucun2008
SSE2High-end Intel200020001999
SSE2AMD20032004 (Barton)2004 (Paris)
SSSE3Atom2008 (first Atom)Aucun2008
SSSE3High-end Intel2006 (Core)2007 (Yonah)2006/2007
SSSE3AMD2011 (Bulldozer, Bobcat)2012 (Sempron X2 198, A4-3450)2012
x86_64Atom200820122013
x86_64High-end Intel2005 (Prescott)2007 (Yonah)2006/2007
x86_64AMD20032004 (Barton)2004
SSE4.2Atom20132013 (Saltwell / Avoton)2013/2014
SSE4.2High-end Intel2008 (Nephalem)2008 (Core)2008
SSE4.2AMD2011 (Bulldozer)2011/2012 (A4-3450)2011/2012
AVXAtomProbablement 2022FuturFutur
AVXHigh-end Intel2011 (Sandy Bridge)2021 (Comet Lake-S)2020/2021 (Tiger Lake)
AVXAMD2011 (Bulldozer)2011/2012 (A4-3450)2011/2012
AVX2AtomProbablement 2022FuturFutur
AVX2High-end Intel20132021 (Comet Lake-S)2020/2021 (Tiger Lake)
AVX2AMD2015 (Carrizo / A10-8700P)2016 (Godavari/A10-7890K)2016
AVX-512AtomAucunAucunAucun
AVX-512High-end Intel2016 (Knights Landing)Futur ou jamaisFutur ou jamais
AVX-512AMDAucunInconnuInconnu
Tableau 1 : support des jeux d’instructions par les processeurs Intel et AMD; High-end Intel correspond à toutes les plateformes d’Intel non basées sur Atom, c’est-à-dire, Xeon, core i3/i5/i7, Pentium Gold, Celeron.

Le tableau ci-dessous est une version simplifiée du précédent, regroupant les trois familles de processeurs : Atom, High-end Intel et AMD.

Jeu d’instructionsDate d’introductionDernier proc à ne pas gérerDébut du support complet sur toute la gammeDélai entre introduction et adoption généralisée
SSE220002004 (Barton)2004 (Paris)4 ans
SSSE32006 (Core)2012 (Sempron X2 198, A4-3450)20126 ans
x86_6420032012201310 ans
SSE4.22008 (Nephalem)2013 (Saltwell / Avoton)2013/20145-6 ans
AVX2011 (Sandy Bridge)FuturFutur> 11 ans
AVX22013FuturFutur> 9 ans
AVX-5122016 (Knights Landing)Futur ou jamaisFutur ou jamais> 6 ans
Tableau 2 : introduction et adoption généralisée des jeux d’instructions sur l’ensemble de la gamme de processeurs Intel et AMD

Rythme d’adoption

On peut observer une adoption très lente de certains jeux, pour diverses raisons. Le jeu d’instructions x86_64 a été introduit en 2003 par AMD sur la gamme Athlon 64 et a été rapidement généralisée par AMD; lorsqu’Intel a réalisé l’échec de l’Itanium, il a adopté cette architecture 64 bits en 2005 avec le Pentium 4 Prescott, ce qui peut paraître assez rapide. Malheureusement, Intel a énormément tardé à passer à l’architecture 64 bits sur la plateforme Atom. L’adoption des autres jeux d’instructions SSE2, SSSE3 et SSE4.2 a été plus rapide, avec un certain retard pour les processeurs Atom ou pour AMD; possiblement à cause d’une certaine concurrence entre le SSE4a d’AMD et les SSSE3/SSE4.1/SSE4.2 d’Intel.

Le cas d’AVX et AVX2 mérite qu’on s’y attarde. AMD a très vite adopté AVX, presque simultanément avec Intel, se basant sur les plans annoncés d’Intel. C’est en se basant sur une ébauche trop précoce qu’AMD s’est retrouvé à implémenter FMA4 alors qu’Intel adopta finalement FMA3. Pourtant, plus de dix ans plus tard, une partie importante des processeurs commercialisés par Intel en 2020 ou 2021 ne gère toujours pas ce jeu d’instructions : Celeron et Pentium Gold de l’architecture Comet Lake ainsi que tous les produits dérivés de la microarchitecture Atom Tremont. Cette situation semble explicable par une fragmentation du marché voulue par Intel, avec une désactivation des unités AVX du Comet Lake sur le bas de gamme (Celeron et Pentium Gold) tandis que les processeurs montant en gamme, basés sur la même microarchitecture (core i3/i5/i7/i9) en bénéficient pleinement.

Point de vue du programmeur

Un développeur d’application peut ignorer les utilisateurs disposant d’une machine très ancienne, car ils représentent une faible fraction de l’ensemble des utilisateurs et pourraient éventuellement avoir des problèmes de performances d’exécution trop sévères pour que l’application soit utilisable. Les hypothèses que l’on peut raisonnablement faire dépendent de la catégorie d’application et de son exigence de performances. Un jeu vidéo dernier cri exigera une machine récente, ce qui est attendu des joueurs. Une application bureautique telle qu’un traitement de texte devra tourner sur de vieilles machines. En 2021, au CHU de Rouen, on trouve encore beaucoup de Pentium G640 datant de 2012, voire plus ancien. Ce retard technologique peut être expliqué par plusieurs facteurs. Premièrement, les constructeurs de PC ont toujours du retard par rapport aux processeurs vendus; pour Lenovo, on peut compter environ un an de retard. Deuxièmement, une grande entreprise préfèrera avoir un parc homogène de machines, quitte à avoir des machines un peu plus anciennes. Troisièmement, le renouvellement du parc est une opération coûteuse qui sera faite d’autant plus rarement que les machines ont une faible obsolescence pour un usage bureautique. Ces machines, veilles de 9 ou 10 ans restent encore suffisamment performantes pour un usage bureautique confortable. En raison de la crise COVID-19 et de l’essor des cryptomonnaies, la demande d’électronique a beaucoup augmenté, les ruptures de stock sont nombreuses et actuellement des machines obsolètes sont vendues hors de prix. C’est ainsi qu’on peut trouver, des machines neuves basées sur des processeurs qui ont 5 ans de retard.

Aujourd’hui on peut raisonnablement considérer que le SSE2 est disponible partout, même pour une application bureautique. On comprend aussi qu’exiger AVX est déraisonnable, même pour un jeu vidéo; ces derniers étant généralement plus exigants sur le GPU que le CPU. C’est ainsi qu’un des jeux les plus gourmands en ressources de l’année 2020, tourne finalement avec un CPU de milieu de gamme de 2012 mais était à deux doigts de ne pas tourner sur l’entrée de gamme de 2020 à cause d’AVX, pour un gain de performances tellement négligeable que les développeurs ont préféré complètement supprimer le support plutôt que de l’utiliser conditionnellement à sa présence. Il faut dire que de gérer optionnellement un jeu d’instructions est très difficile à faire correctement à cause des fonctions inlines, sauf sur le compilateur C d’Intel (ICC) qui pénalise abusivement les processeurs AMD en exécutant le pire code dessus. Pour compliquer les choses, l’AVX ne fait pas bon ménage avec le SSE/SSE2, parce que les instructions SSE opérant sur les registres xmm préservent la partie haute des registres ymm; cette fonctionnalité étrange étant voulue par Intel pour des raisons de rétro-compatibilité avec des pilotes binaires susceptibles d’exécuter des instructions SSE durant une interruption matérielle. Cela engendre une pénalité énorme si VZEROUPPER n’est pas appelé entre les instructions AVX et SSE; sauf pour les Xeon Phi de la microarchitecture Knights Landing pour lesquels c’est exactement l’inverse, avec une pénalité de 36 cycles à exécuter VZEROUPPER. On peut raisonnablement ignorer les Xeon Phi qui concernent une toute petite partie du marché, mais l’obligation d’employer VZEROUPPER est déjà délicate dans tout programme utilisant des fonctions de bibliothèques basées sur SSE/SSE2 au milieu de code AVX. Ce problème serait transitoire si le support d’AVX avait été rapidement universel, mais Intel n’a pas voulu que ce soit ainsi. Heureusement, mes analyses semblent montrer que même pour des opérations très parallélisables, AVX+AVX2+FMA3 n’apportent quasiment aucune accélération par rapport à SSE2 : environ 5% pour la multiplication de matrices sur Ryzen 1700. En bref, AVX et AVX2 sont des jeux d’instructions utiles pour le marketing mais trop compliquées à utiliser, trop peu disponibles et ne servant finalement pas à grand chose. C’est encore pire avec AVX-512 qui finalement consomme tellement d’espace et d’énergie qu’il ne sera pas disponible sur la prochaine génération de processeurs Alder Lake d’Intel parce que ces derniers utilisent une stratégie big.LITTLE.

Robuste ou pas ?

Vous connaissez peut-être le test de Brown-Forsythe publié par Brown et Forsythe et décrit par les auteurs comme un test d’égalité des variances robuste aux écarts à la normalité. Le test est très simple à comprendre : il s’agit d’une ANOVA comparant la moyenne des écarts à la médiane. S’il y a seulement deux groupes, c’est équivalent à un test de Student sur la variable z égale aux écarts à la médiane du groupe.

Il y a un problème avec la manière de présenter les choses. De deux choses l’une:

  1. Soit, c’est un test de comparaison des écarts moyens à la médiane. Auquel cas, il est assez robuste aux écarts à la normalité.
  2. Soit, c’est un test de comparaison de variances, extrêmement fragile, reposant sur certaines hypothèses douteuses et dont le risque alpha peut facilement frôler les 100%

Considérons le cas de figure suivant: considérons comme première distribution une loi exponentielle de paramètre lambda=1. Cette loi exponentielle a une variance et une moyenne toutes deux égales à 1. Considérons comme seconde distribution une loi normale centrée-réduite de moyenne 0 et de variance 1. Ces deux distributions ont exactement la même variance et correspondent donc à l’hypothèse nulle d’égalité des variances du test. Pourtant, dès que les échantillons sont de taille non négligeable, le test va rejeter cette hypothèse nulle avec une probabilité très forte:

Le code R ci-dessous illustre le problème :

library(lawstat)
set.seed(2021)
pvalues=replicate(1000, {
	N=1e4
	y=c(rnorm(N), rexp(N))
	group=rep(c(0,1), c(N, N))
	lawstat::levene.test(y, group, location="median")$p.value
})
max(pvalues) # renvoie 8e-14

Ainsi, bien qu’on soit sous l’hypothèse nulle, sur 1000 tests, on observe 1000 tests significatifs, soit un risque alpha supérieur à 99%. Pourquoi ce résultat aussi lamentable ? Parce que le test de Brown-Forsythe rejette en réalité l’égalité de l’écart moyen à la médiane : soit environ 0,80 pour la loi normale centrée-réduite et 0,69 pour la loi epxonentielle. Dans la même idée, on peut aussi avoir une puissance à 5% si les variances sont inégales mais que les écarts moyens à la médiane sont égaux.

Le test de Brown-Forsythe est un test de comparaison de variances valide sous une hypothèse très forte, c’est-à-dire, une hypothèse qui, si elle n’est pas vérifiée, rend très probable une inflation du risque alpha frôlant les 100% : la forme des distributions doit être totalement identique et seule la variance doit différer entre les deux distributions. C’est-à-dire, que si X et Y sont deux variables aléatoires suivant respectivement les deux lois dont ont compare les variances, et s’il existe deux constantes a et b tels que distribution(X)=distribution(a*Y+b), alors le test de Brown-Forsythe contrôle le risque alpha. On repose donc sur une hypothèse « shift in location and scale ». Cette hypothèse est suffisante à la validité du Brown-Forsythe, mais au sens strict, n’est pas tout à fait nécessaire, puisqu’il suffit que la famille de distributions soit choisie telle que l’égalité des écarts moyens soit équivalente à l’égalité des variances. En pratique, l’hypothèse « shift in location and scale » est une bonne approximation de la condition de validité du Brown-Forsythe, qui devient alors un test extrêmement fragile.

On peut argumenter que ce n’est pas grave : le test compare les « dispersions » et pas les variances au sens strict. Seulement, ce n’est pas présenté comme ça par les auteurs et la distinction a de l’importance. Par exemple, le test de Student repose sur l’hypothèse d’égalité des variances, pas d’égalité des écarts moyens à la médiane. La vérification de l’hypothèse d’homoscédasticité du test de Student par un test d’égalité des variances pose plus de problèmes qu’elle n’en résout, mais c’est une autre histoire.

Mann-Whitney

Dans le même ordre d’idée, le test de Mann-Whitney est souvent interprété, à tort, comme un test de comparaison de médianes. Il n’est valide, interprété ainsi, qu’en condition d’hypothèse de shift-in-location, c’est-à-dire qu’une distribution doit être simplement décalée d’une constante par rapport à l’autre : distribution(X)=distribution(a+Y). Sur un gros échantillon, le test de Mann-Whitney tendra vers un risque alpha à 100% pour la comparaison de la distribution exponentielle avec lambda=log(2) dont la médiane est à 1, et la loi normale de moyenne (et médiane) égale à 1 et de variance égale à 1. Sous l’hypothèse « shift-in-location », l’égalité des moyennes est strictement équivalente à l’égalité des médianes, et le test de Mann-Whitney est autant un test d’égalité des médianes qu’un test d’égalité des moyennes. Ci-dessous l’exemple:

wilcox.test(rexp(1e4,rate=log(2)), 1+rnorm(1e4))

Wilcoxon sur séries appariées

Le test de Wilcoxon sur séries appariées n’est pas non plus un vrai test d’égalité des médianes, sauf sous l’hypothèse shift-in-location : il s’agirait plutôt d’un test de comparaison de la pseudo-médiane des différences appariées à la valeur zéro. Ci-dessous l’exemple illustrant le propos, avec un risque alpha proche de 100%:

wilcox.test(rexp(1e4,rate=log(2)),  1+rnorm(1e4), paired=TRUE)

Théorisation de la problématique

On peut grossièrement catégoriser les biais des tests statistiques et estimateurs statistiques en trois catégories, en cas d’écart à leurs hypothèses :

  1. Ceux qui tendent à diminuer avec la taille d’échantillon et finissent par s’annuler (asymptotiquement non biaisé)
  2. Ceux qui sont indépendants de la taille d’échantillon (ou presque)
  3. Ceux qui tendent à croitre avec la taille d’échantillon pour se rapprocher du biais maximal théoriquement imaginable (p.e. risque alpha à 100%)

En théorie, toute relation est possible entre les biais et la taille d’échantillon, mais en pratique, vous devriez presque toujours rencontrer des biais rentrant dans ces catégories si vous utilisez des outils de conception pas trop tirée par les cheveux.

Les biais engendrés par l’hétéroscédasticité pour le test de Student sur des échantillons de taille inégale, rentrent dans la catégorie N°2. Le biais engendré par les écarts à la normalité de la distribution, mais avec homoscédasticité, pour le test de Student, rentre dans la catégorie N°1. Comme précédemment présenté, en cas d’écarts aux hypothèses shift-in-location du test de Wilcoxon sur séries appariées ou Mann-Whitney, lorsque ces tests sont interprétés comme des tests de comparaison de médianes, les biais engendrés rentreront généralement dans la catégorie N°3.

Étant donné que toutes les hypothèses faites en biostatistiques sont plus ou moins fausses, je n’accepte jamais les biais de la catégorie N°3. C’est-à-dire, qu’il vaut mieux changer d’estimateur que d’accepter l’inacceptable. Avec certaines méthodes très génériques comme le boostrap, on peut totalement éliminer les biais des catégories N°1 et N°2 de toute façon.

Les biais de la catégorie N°1 peuvent être préoccupants sur de petits échantillons mais rarement sur de grands échantillons. Par exemple, pour le test de Student sur une variable écartée de la loi normale mais restant dans un intervalle borné, il n’y a pas de souci à se faire dès qu’on dépasse les quelques centaines d’observations dans chaque groupe. Enfin, les biais de la catégorie N°2 doivent faire l’objet d’une grande prudence. Ils peuvent être acceptables lorsqu’on sait que les déviations aux hypothèses sont modestes, mais si on a un doute, mieux vaut employer un estimateur asymptotiquement correct comme le boostrap, c’est-à-dire, dont tous les biais vont tendre vers zéro sur de grands échantillons.

C’est une des raisons pour lesquelles je suis friand de bootstrap, avec ses autres propriétés séduisantes : le boostrap est une méthode très générale d’estimation, fonctionnant sur toute statistique customisée, reposant sur très peu d’hypothèses (juste des observations indépendantes et identiquement distribuées), présentant des biais asymptotiquement nuls, sauf exception, et présentant des outils de diagnostics pertinents. Par ailleurs, il oblige à conceptualiser la méthodologie d’échantillonnage qui est souvent ignorée. Bref, le boostrap, mangez-en ! C’est bon ! Enfin, faut-il encore le faire correctement, parce faire du boostrap avec de la stepwise regression mais appliquer le bootsrap après la sélection des variables, c’est du grand n’importe-quoi mais ça se voit souvent.

Exemple d’optimisation du code en R

Un billet pour expliquer les étapes conduisant à une amélioration des performances d’un programme assez simple : réalisant un bootstrap résiduel dans un modèle linéaire général. Cette méthode statistique est adaptée aux écarts assez sévères à la normalité, fragilisant le théorème central limite et limitant la validité de la méthode de Wald. Cette méthode ne résout en rien les défauts d’additivité ou de linéarité du modèle ou les problèmes de corrélation entre observations. Dans le scenario à deux groupes, la méthode est équivalente au parabootstrap avec variances égales. On peut la résumer à la procédure suivante pour un échantillon de taille N:

  1. Tirer au hasard N valeurs parmi les N résidus, avec remise.
  2. Ajouter ces N nouveaux résidus aux N valeurs ajustées (fitted values), créant ainsi un nouveau vecteur réponse
  3. Réajuster (refitting) le modèle à ce nouveau vecteur réponse, en conservant la matrice de modèle (model matrix) inchangée
  4. Répéter les étapes (1) à (3) un grand nombre de fois (p.e. dix mille fois) tout en enregistrant les coefficients du modèle à chaque refitting
  5. Calculer les intervalles de confiance des coefficients comme les quantiles des coefficients de l’ensemble des expériences enregistrées en (4)

L’échantillon de test était issu d’une vraie étude, et comportait 924 observations avec 11 covariables qualitatives conduisant à 38 coefficients distincts après recodage des variables qualitatives à k modalités en k-1 variables binaires et ajout d’une ordonnée à l’origine (intercept).

La première approche est directe, basée sur le package ‘boot’ et sur la fonction ‘lm’:

resid_boot1=function(mod, R1=10, R2=1000) {
	f=fitted(mod)
	res=resid(mod)
	boot(data=model.frame(mod), statistic=function(data, idx) {
		datax=data;
		# on change la variable réponse
		datax[,1]=f+res[idx]
		# puis on réajuste (refitting)
		coef(lm(data=datax, formula=formula(mod)))
	}, R=R1*R2)
}

Le paramètre mod représente un modèle ajusté (fitted) sur le jeu de données original et R1 et R2 sont deux paramètres représentant le nombre de simulations. Habituellement, on utilise un seul paramètre R représentant le nombre total de simulations, mais pour avoir une homogénéité d’interface, j’ai ici utilisé deux paramètres dont le produit est égal au nombre total de simulations. On pourra remarquer que pour réajuster (refitting) le modèle, je n’ai pas fait appel à la fonction update() de R. Cette fonction est extrêmement dangereuse et ne devrait jamais être utilisée parce qu’elle ne réajuste pas le modèle dans l’environnement dans lequel il avait été ajusté la première fois.

Une fois cette fonction exécutée, il suffit d’appeler boot.ci() pour chacun des coefficients du modèle afin de calculer les intervalles de confiance:

bt = resid_boot1(original_model, R1=10,R2=100)

confidence_intervals = t(sapply(seq_along(bt$t0),
function (i) {
  boot.ci(bt, type="perc",index=i)$perc[4:5]
}))

La fonction resid_boot1 est malheureusement lente pour R version 4.0.3. Sur un PC assez moderne (AMD Ryzen 7 1700, 3 Ghz, 8 coeurs, 24 Go de RAM, bien ventilé, tournant en turbo à 3.3 Ghz en continu) sous Windows 10 64 bits il faut 4,8 secondes pour mille simulations (R1=10, R2=100).

Le premier constat qui permet d’améliorer les performances, c’est que l’appel à la fonction lm() refait itérativement des calculs inutiles. Cette fonction va commencer par recoder les variables explicatives qualitatives en variables binaires, mais comme ces variables explicatives ne changent pas d’une simulation à l’autre, le recodage est redondant. Pour gagner en temps, on peut utiliser la fonction lm.fit() qui prend directement en entrée, une matrice numérique dite matrice de modèle (model matrix) qui a déjà toutes les colonnes recodées.

Ainsi, la deuxième implémentation du boostrap résiduel se base sur lm.fit:

resid_boot2=function(mod, R1=10, R2=1000) {
	f=fitted(mod)
	res=resid(mod)
	boot(data=model.matrix(mod), statistic=function(data, idx) {
		lm.fit(data, f+res[idx])$coefficients
	}, R=R1*R2)
}
system.time(bt2<-resid_boot2(modT, R1=10, R2=100))

Pour le même nombre de simulations (mille), il ne faut plus que 1,42 secondes soit une multiplication par 3,4 de la vitesse d’exécution.

Cette fonction lm.fit fait quelques vérifications avant d’appeler .lm.fit, qui fonctionne presque de la même manière mais devrait être un peu plus rapide; avec cette dernière fonction le temps d’exécution passe à 1,36 seconde; amélioration minime.

L’implémentation de .lm.fit de R n’étant pas optimale, il existe une fonction lmfit dans le package Rfast, dont l’interface est semblable mais qui devrait être un peu plus rapide. Le temps d’exécution passe à 0,97 seconde avec Rfast::lmfit.

L’étape suivante d’optimisation repose sur une propriété importante du modèle linéaire général estimé par les moindres carrés: conditionnellement à la matrice de modèle (c’est-à-dire aux variables explicatives), chacun des coefficients du modèle n’est qu’une combinaison linéaire (c’est-à-dire, une moyenne pondérée) du vecteur réponse. Cela se voit très bien dans la formule classique d’estimation du modèle linéaire général. Étant donné une matrice de modèle X à n lignes (représentant les n observations) et k colonnes (représentant les k variables explicatives après recodage) et un vecteur réponse Y à n valeurs (la variable à expliquer), les coefficients du modèle s’estiment comme (X’ X)-1X’ Y où X’ représente la matrice transposée de X. Ainsi, on commence par faire le produit de la transposée de X avec X elle même, puis on en prend l’inverse; ensuite on remultiplie par la transposée de X; enfin, on multiple par le vecteur Y. Seule cette dernière étape va varier dans le boostrap résiduel, puisque X est constante. On peut donc gagner une inversion de matrice et trois produits de matrice en précalculant (X’ X)-1X’.

C’est donc l’étape d’optimisation suivante:

resid_boot5=function(mod, R1=10, R2=1000) {
	f=fitted(mod)
	res=resid(mod)
	x = model.matrix(mod)
	linmat = solve(t(x) %*% x) %*% t(x)
	boot(data=model.matrix(mod), statistic=function(data, idx) {
		(linmat %*% (f+res[idx]))[,1]
	}, R=R1*R2)
}

Le temps d’exécution passe à 0,12 seconde pour mille itérations, soit environ huit fois plus rapide que la version précédente; afin de pouvoir calculer ce chiffre avec précision, il faut multiplier par dix le nombre de simulations : 1,20 seconde pour dix mille simulations.

Les optimisations suivantes sont des « micro optimisations » qui feront quand même gagner une vistesse non négligeable (vous verrez). Afin de mieux calculer les timings et de se défaire de quelques propriétés inimportunes du package boot sur la manière de générer les nombres aléatoires (en un énorme bloc si simple=FALSE ou en trop petits blocs si simple=TRUE), on peut réécrire la même fonction de manière plus brute. Au passage, on peut transposer la matrice avant la multiplication de matrice et utiliser crossprod qui est un peu plus rapide. En même temps, on va réajuster (refitter) R1=10 simulations à chaque itération, afin de réduire les coûts intertiels (overheads) du produit matriciel et du tirage au sort des observations.

resid_boot6=function(mod, R1=10, R2=1000) {
	f=fitted(mod)
	res=resid(mod)
	ff = matrix(fitted(mod), nrow=nrow(x), ncol=R1)
	x = model.matrix(mod)
	tlinmat = t(solve(t(x) %*% x) %*% t(x))
	N=nrow(x)
	a = replicate(R2, crossprod(tlinmat,
                                    ff+sample(res,N*R1,replace=TRUE)))
	dim(a)=c(length(coef(mod)),R1*R2)
	list(t=t(as.matrix(a)), t0=coef(mod))
}

Pour dix mille itérations, on passe à 0,97 seconde, ce qui est un maigre gain par rapport aux 1,20 secondes de l’implémentation précédente.

L’étape suivante a consisté à évaluer si le goulot d’étranglement (bottleneck) se situait dans le produit matriciel ou l’échantillonnage des résidus, en décomposant la fonction. On pouvait constater que l’échantillonnage des résidus était le facteur limitant, comme si le générateur pseudo-aléatoire était trop lent. En essayant tous les générateurs pseudo-aléatoires proposés par R (cf fonction set.seed), j’observais une vitesse d’exécution presque identique, mais lente, pour tous les générateurs. Même en se basant sur la fonction runif(), qui est la plus directe des fonctions, on ne dépassait pas les 15,3 millions de nombres aléatoires générés par seconde, comme le montre la ligne de code suivante:

system.time(for(i in 1:1e4) sample.int(1:nobs(model),1e4,replace=TRUE))
#   user  system elapsed 
#   6.55    0.00    6.55

Pour un processeur qui tourne en permanence en mode turbo à 3,3 Ghz, ça fait environ 216 cycles de CPU pour calculer un nombre aléatoire. Pour l’algorithme Mersenne-Twister par défaut, ça paraît beaucoup, puisque cet algorithme est connu pour être plutôt rapide, basé sur quelques xor et bit shiftings par nombre 32 bits généré. Afin de vérifier la vitesse optimale de génération des nombres aléatoires selon cet algorithme, j’ai téléchargé une version très basique de l’algorithme (https://github.com/ESultanik/mtwister) et l’ai compilée avec GCC 8.3.0 pour x86_64 et l’option d’optimisation -O2. Elle générait un milliard de nombres pour 3.9 seconde, soit à peine 13 cycles de CPU par itération. Cela correspondait globalement à une multiplication par 16,7 de la vitesse de génération des nombres aléatoires. C’est ainsi qu’il devenait rentable de réimplementer une fonction d’échantillonnage (fastsample) en C:

Le cœur de l’algorithme est donné ci-dessous en code C:

	for(R_xlen_t i=0; i<inobs; i++) {
		double unif=(((double)genRandLong(rangen))
                            * (1.0/(1+(double)0xFFFFFFFFL)));
		presult[i] = psamp[(size_t)(ssize*unif)];
	}

On pouvait alors remplacer sample() par fastsample() pour quasiment annihiler le temps de génération des nombres pseudo-aléatoires:

resid_boot7=function(mod, R1=10, R2=1000, seed=1337) {
	rgen=rangen(seed)
	f=fitted(mod)
	res=resid(mod)
	ff = matrix(fitted(mod), nrow=nrow(x), ncol=R1)
	x = model.matrix(mod)
	tlinmat = t(solve(t(x) %*% x) %*% t(x))
	N=nrow(x)
	a = replicate(R2, crossprod(tlinmat, ff+fastsample(rgen, res,N*R1)))
	dim(a)=c(length(coef(mod)),R1*R2)
	list(t=t(as.matrix(a)), t0=coef(mod))
}

Pour dix mille itérations, il ne fallait plus que 0,40 seconde. Sur ces 0,40 seconde, le fastsample représentait environ 0,04 seconde. Pour ce calcul, j’ai dû multiplier par dix le nombre d’itérations (soit cent mille itérations). Presque tout le temps de calcul restant provenait de la multiplication de matrice (crossprod). Cette multiplication, bien que censée être optimisée, est loin d’être optimale comme on peut le voir dans la fonction internal_crossprod du fichier src/main/array.c du code source de R. C’est une implémentation assez naïve, basée sur les bonnes vieilles FP80 du FPU x87. En utilisant les instructions vectorielles (SIMD) avec des FP64, on peut a priori gagner en performance. Pour bien utiliser les SIMD FP64 des processeurs x86, c’est-à-dire, les SSE2 et AVX/AVX2, il faut aligner les accès à la mémoire sur 128 bits (SSE2/AVX) ou 256 bits (AVX2). Cela oblige à arranger les matrices pour ces alignements. Par exemple, une matrice 7×3 (7 lignes et 3 colonnes) va devoir être transformée en une matrice 8×3, avec la dernière ligne à zéro pour que les accès mémoire soient bien alignés, sans compter l’allocation du vecteur de données qui devra être bien alignée. Cela va demander donc un petit temps de transformation de la matrice; heureusement, pour la matrice (X’X)-1X’, ce travail peut n’être fait qu’une seule fois. C’est ainsi, qu’il a été possible d’écrire deux fonctions en langage C en utilisant les Vector Extensions de GCC: une fonction plmat() qui prépare la matrice à gauche de la multiplication (transposition/alignement) et une fonction fastmatmul() qui multiplie cette matrice préparée par une nouvelle matrice. En utilisant les SSE2, la nouvelle fonction s’exécutait en 0,14 seconde pour dix mille itérations. En passant à cent mille itérations (R1=100, R2=1000), il fallait 1,24 seconde. Avec AVX2+FMA3, on passait à 1,18 seconde pour cent mille itérations. À ce stade, il ne sert à rien d’optimiser par rapport à un usage ponctuel du bootstrap résiduel. Cependant pour évaluer les propriétés statistiques (notamment les biais) du boostrap résiduel, on peut avoir besoin de faire des dizaines de milliers de bootstrap. En visant dix mille simulations par bootstrap, on peut facilement se retrouver avec des centaines de millions de calculs. C’est là qu’on peut gagner en performance en utilisant les multiples cœurs du processeur; cela se fait bien avec le package ‘snow’ sous Windows ou parallel::mclapply sous GNU/Linux. En pratique, avec huit cœurs on gagne sans difficulté un facteur huit; avec le SMT, on peut encore multiplier la vitesse par 1,15 ou 1,20 environ, mais au-delà il est difficile de gagner en performances sur un CPU. Sur GPGPU on peut espérer y gagner encore. Cependant mes résultats préliminaires avec OpenCL sur une GeForce GTX 1060 3Go sont extrêmement décevants, avec des performances pires que sur CPU avec OpenCL: cela semble s’expliquer par un usage trop intensif de nombres entiers. Il y a probablement des optimisations auxquelles je n’ai pas pensé et qui pourraient accélérer.

Le tableau ci-dessous récapitule les gains de performance successifs:

AlgorithmeTemps estimé pour 100 000 itérations (secondes)
1 : naïf avec lm()480
2 : amélioré avec lm.fit()142
3 : amélioré avec .lm.fit()136
4 : avec Rfast::lmfit97
5 : précalculant (X’X)-1X’11,9
6 : pré-transposition et calcul de 10 itérations simultanées9,7
7 : réécriture fastsample() en C4,0
8 : optimisation multiplication de matrices avec SSE20,124
9 : optimisation AVX2+FMA30,118
10 : OpenCL+Geforce GTX 1060 3 Go0,13
Performance des différents algorithmes de bootstrap résiduel du modèle linéaire général

La version N°9 est environ quatre mille fois plus rapide que la version N°1. Il n’y avait qu’une seule vraie optimisation algorithmique (méthode N°5), faisant gagner un facteur 8, les autres optimisations étant des petites bidouilles qui ne changent pas la complexité algorithmique.

Pseudo-aléatoire R

Pourquoi le générateur pseudo-aléatoire est-il aussi lent ? Même l’accès le plus direct à la source pseudo-aléatoire, la fonction runif() est trois fois plus lente que fastsample() bien qu’elle ne fasse pas l’indexation d’un vecteur. En fouillant un peu le code source de R, en appelant les différentes parties pour en évaluer les performances et en regardant le code assembleur généré par GCC, on découvre que le problème ne vient pas du générateur Mersenne-Twister lui-même (fonction MT_genrand dans src/main/RNG.c) ni de unif_rand dans src/main/RNG.c dont les inerties (overheads) sont faibles. Le problème vient de random2() dans src/library/stats/src/random.c. Soit dit au passage, le code de src/main/random.c est redondant avec celui de src/library/stats/src/random.c. Le premier code ne semble finalement jamais être appelé en situation ordinaire et il est plus performant que le second.

> runif
function (n, min = 0, max = 1) 
.Call(C_runif, n, min, max)
<bytecode: 0x0000000004706130>
<environment: namespace:stats>
> system.time(for(i in 1:1e4)runif(1e4))
   user  system elapsed 
   1.48    0.00    1.49
> system.time(for(i in 1:1e4).Internal(runif(1e4,0,1)))
   user  system elapsed 
   1.28    0.00    1.29

Sur l’AMD Ryzen 1700, on constate donc des performances un peu meilleures pour la fonction interne, implémentée dans src/main/random.c que la fonction de bibliothèque dynamique implémentée dans src/library/stats/src/random.c, même si les deux restent mauvaises. Concernant src/library/stats/src/random.c on comprend rapidement le problème en lisant le code source:

for (R_xlen_t i = i0; i < n; i++) {
    rx[i] = fn(ra[i % na], rb[i % nb]);
    if (ISNAN(rx[i])) naflag = TRUE;
}

Où la fonction fn() génère un nombre aléatoire selon une loi uniforme; les deux paramètres de fn représentent respectivement le minimum et le maximum de la loi uniforme.

Le problème de performance provient des deux divisions entières (i % na et i % nb) réalisées pour le cyclage des arguments min et max de runif pour gérer des appels tels que stats::runif(100, 1:5, c(50,100)). La division entière est une opération très lente sur de nombreux processeurs, nécessitant 14 à 45 cycles (reciprocal throughput) sur la microarchitecture Zen 1 (https://www.agner.org/optimize/instruction_tables.pdf). Sur processeur Intel, c’est encore pire; sur Coffe Lake (8ème génération d’Intel Core) la latence est de 42 à 95 cycles et le reciprocal throughput est de 24 à 90 cycles. Le problème est moins sévère avec le code de src/random qui utilise MOD_ITERATE2_CORE (src/include/R_ext/Itermacros.h), dont le code est montré ci-dessous:

define MOD_ITERATE2_CORE(n, n1, n2, i, i1, i2, loop_body) do { \
        for (; i < n;                                                   \
             i1 = (++i1 == n1) ? 0 : i1,                                \
                 i2 = (++i2 == n2) ? 0 : i2,                            \
                 ++i) {                                                 \
            loop_body                                                   \
                }                                                       \
    } while (0)

Le code assembleur généré utilise quand même des cmov suboptimaux et a du mal à tenir dans les registres. On peut accélérer énormément les choses en réécrivant le code pour le cas particulier où les arguments min et max à runif() sont des vecteurs de longueur 1. Le nouveau code de random2() dans src/main/random.c devient alors:

static Rboolean random2(double (*f) (double, double),
                        double *a, R_xlen_t na, double *b, R_xlen_t nb,
                        double *x, R_xlen_t n)
{
    double ai, bi;
    R_xlen_t i, ia, ib;
    Rboolean naflag = FALSE;
    errno = 0;
    if (na == 1 && nb == 1) {
        double ra=*a, rb=*b;
        double *end = x+n;
        double *x0=x;
        while(x < end) {
            *x = f(ra, rb);
            if (ISNAN(*x)) naflag = TRUE;
            x++;
        }
        *x0=(R_xlen_t)f;
    } else {
        MOD_ITERATE2(n, na, nb, i, ia, ib, {
            ai = a[ia];
            bi = b[ib];
            x[i] = f(ai, bi);
            if (ISNAN(x[i])) naflag = TRUE;
        });
    }
    return(naflag);
}

Par ailleurs la fonction runif dans src/nmath/runif.c qui est chargée de calculer un seul nombre aléatoire entre min et max, est sous-optimale dans le cas où min=0 et max=1, parce qu’elle va faire une addition, une soustraction et une multiplication inutiles pour fournir un résultat compris entre min et max. On peut alors gagner un petit peu à tester le scenario commun min=0 et max=1:

double runif(double a, double b)
{
    if (a == 0 && b == 1) { /* ajout pour accélérer */
        return unif_rand();
    }
    if (!R_FINITE(a) || !R_FINITE(b) || b < a)  ML_WARN_return_NAN;

    if (a == b)
        return a;
    else {
        double u;
        /* This is true of all builtin generators, but protect against
           user-supplied ones */
        do {u = unif_rand();} while (u <= 0 || u >= 1);
        return a + (b - a) * u;
    }
}

Moyennant ces modifications, les appels à runif() bénficièrent d’une accélération par un facteur trois environ. Il existe encore quatre ralentissements non négligeables que je pus évaluer en réécrivant complètement la fonction pour shunter certaines étapes : (1) l’usage d’une boucle générique pour tous les générateurs pseudo-aléatoires à deux paramètres (rnorm, rbinom, runif, etc.) engendre un niveau d’indirection (2) le test ISNAN(*x) répétitif n’est pas nécessaire pour une distribution uniforme et engendre un UCOMISD et un CMOVP à chaque itération (3) la comparaison a==0 && b == 1 dans runif(double a,double b) est exécutée à chaque itération (4) la fonction unif_rand() qui s’occupe de sélectionner le générateur pseudo-aléatoire ajoute encore un niveau d’indirection par rapport à un appel direct à MT_genrand() pour l’usage du générateur Mersenne-Twister. Cette dernière sélection du générateur est particulièrement peu performante, car basée sur un gros switch() après lecture d’une variable globale. Si on corrige ces 4 problèmes on gagne encore un facteur deux en performances. C’est ainsi qu’on peut multiplier par six les performances de runif.

En mettant ça au propre, je pense que je pourrais écrire un petit patch pour R.

Tests anti-hiérarchiques

Vous connaissez peut-être les tests hiérarchiques; une manière de gérer la multiplicité des tests en planifiant une liste de tests, à réaliser dans un ordre précis. Chaque test est réalisé à un niveau de significativité donné (5% le plus souvent). Les tests sont effectués tant qu’ils sont tous significatifs (p < 5%) et il y a un arrêt de la procédure au premier test échoué (p > 5%).Cela permet de contrôler le risque alpha global à 5% ainsi que les risques alpha associés à chaque test.

Voici la procédure en résumé:

  1. Une liste de tests est pré-planifiée
  2. On effectue tous les tests dans la liste, tant que tous les résultats sont significatifs (p < 5%)
  3. On s’arrête de tester au premier résultat non significatif (p > 5%)

J’ai constaté que dans la pratique, une approche inversée est fréquemment employée:

  1. La liste de tests n’est pas pré-planifiée, mais est décidée au fur et à mesure, selon les résultats précédents
  2. On effectue des tests tant que les résultats sont tous non-significatifs (p > 5%)
  3. On s’arrête de tester au premier résultat significatif (p < 5%)

Vous comprenez alors pourquoi j’appelle cela les tests anti-hiérarchiques.

Cette méthode est capable d’engendrer une inflation majeure du risque alpha malgré un Bonferroni, parce que le nombre de tests n’est pas défini à l’avance. En effet, même si on applique rigoureusement un Bonferroni à chaque étape, le risque alpha sera environ égal à 5% (1er test) + 5%/2 (2ème test) + 5%/3 (3ème test) + 5%/4 + …

En réalité, cette formule est inexacte car les risques ne s’additionnent pas tout à fait. En considérant que les risques sont indépendants, on peut calculer qu’avec 20 tests, le risque alpha atteint 16,6% et avec 100 tests, il atteint 23% (cf code R qui suit).

1-prod(1 - 0.05/(1:20))
1-prod(1 - 0.05/(1:100))

Si les tests sont positivement corrélés, le risque croit moins vite que ça. Néanmoins, on ne contrôle toujours pas le risque alpha, malgré un Bonferroni appliqué à chaque nouveau test réalisé.

Il existe pourtant une méthode de contrôle du risque alpha des tests anti-hiérarchiques : réaliser chaque nouveau test avec un niveau de significativité deux fois inférieur au précédent, tout en commençant le 1er test à la moitié du risque alpha.

Ainsi, pour contrôler un risque alpha à 5%, le premier test sera réalisé au seuil de significativité 5%/2, le second sera réalisé au seuil 5%/4, le troisième au seuil 5%/8, le quatrième au seuil 5%/16 et le n-ième au seuil 5%/(2^n). La série mathématique additionnant tous ces seuils tend vers 5%, et donc, quel que soit le nombre de tests réalisés, le risque alpha reste contrôlé.

J’attends avec impatience le 1er article qui contiendra dans sa partie statistique : « nous avons post-planifié l’ensemble des analyses par une méthode anti-hiérarchique avec procédure de correction de multiplicité de type Bonferroni d’acier upgradé »

Analyses en sous-groupes

Un billet que j’écris suite à une remarque d’un reviewer qui a remis en cause ma manière de faire des analyses en sous-groupes, proposant une autre méthode, plus classique mais dont la validité est très discutable. Mes réflexions rejoignent une citation du Pr Bruno FALISSARD que je paraphrase ici, parce que je ne l’ai pas mémorisée de manière exacte : tout le problème est de savoir à qui est la charge de la preuve ?

Ce billet de blog contient du code R. Il n’est pas nécessaire à la compréhension du propos, mais peut vous aider à reproduire les résultats discutés et à comprendre comment ils ont été calculés.

Problématique

Une fois qu’on a mis en évidence le bénéfice d’un traitement sur un grand échantillon de sujets plus ou moins homogènes, cette homogénéité dépendant en partie des critères d’inclusion, on peut raisonnablement craindre que ce bénéfice ne soit pas identique chez tous les individus. Par exemple, on pourrait craindre qu’un traitement bénéfique sur une forme sévère de la maladie, soit inutile voire délétère sur une forme légère, comme on peut l’imaginer pour la dexaméthasone dans le traitement de la COVID-19. On parle alors d’interaction statistique. On distingue les interactions quantitatives des intraction pour lesquelles le traitement a un effet d’amplitude variable selon le sous-groupe mais reste toujours de même signe, des interactions qualitatives pour lesquelles le traitement a des effets opposés selon le sous-groupe: il serait nocif pour certains individus et bénéfique pour d’autres. Ce sont ces dernières interactions qui sont les plus préoccupantes, parce qu’elles remettent en cause la prescription à des sous-groupes identifiables, alors que des interactions quantitatives modestes ne changent pas grand chose aux indications et contre-indications.

Approche N°1 : tests en sous-groupes

La première approche, intuitive, à cette problématique, c’est de répéter le test d’efficacité dans chacun des sous-groupes (p.e. patients âgés, puis patients jeunes). Dans les sous-groupes pour lesquels on arrive à prouver l’effet bénéfique du traitement (p<0.05, IC95% effet = +10% à +30%), on est rassuré et on peut conclure sans crainte à son bénéfice. Dans les sous-groupes pour lesquels on échoue à prouver l’effet bénéfique ou nocif du traitement (p>0.05, IC95% effet=-20% à +30%), on admet qu’on ne peut pas conclure de manière certaine. Dans les sous-groupes pour lesquels on arrive à prouver un effet nocif (p<0.05, IC95% effet = -30% à -10%), on conclut qu’il est contre-indiqué.

Cette approche souffre d’une principale limitation. On manque souvent de puissance pour ces analyses en sous-groupes, conduisant à de nombreux tests non concluants (p>0.05). Si on avait une puissance à 90% avec risque alpha bilatéral à 5% (ou 2,5% unilatéral) pour l’analyse principale, incluant tous les sujets, et qu’on fait une analyse dans 3 sous-groupes parfaitement équilibrés (33% – 33% – 33% de l’échantillon complet), alors, sans même correction de multiplicité des tests la puissance descend à 46% dans chacun des deux sous-groupes, selon la formule R suivante:

pnorm((qnorm(0.975)+qnorm(0.90))/sqrt(3)-qnorm(0.975))

Certains sous-groupes sont souvents plus petits qu’un tiers de l’échantillon complet, et la puissance de l’analyse principale est souvent inférieure à 90%, conduisant à beaucoup de résultats non concluants.

À ce problème, je répondrai que si on n’a pas les données, on ne va pas les inventer. Si on craint réellement à la nocivité d’un traitement dans un tout petit sous-groupe, mais qu’il est trop petit pour qu’on puisse retrouver l’effet dedans, alors on n’a vraiment pas la possibilité de se rassurer, sans lancer d’autres étude sur le sujet.

Variantes

Pour comprendre cette sous-section, il est préférable de lire les chapitres suivants. Revenez-y à la fin de votre lecture.

Doit-on faire une correction de multiplicité des tests ? Ça dépend un peu de l’intention. Si l’analyse principale démarre sur une recherche d’effets dans dix sous-groupes, avec une conclusion portant sur la ou les analyses qui ressortent statistiquement, alors le problème de multiplicité des tests s’applique et une correction paraît pertinente. Il n’est pas nécessaire de corriger la multiplicité des tests si l’analyse principale concerne l’effet global poolé, assumant la « dictature de la moyenne », puis que des analyses en sous-groupes sont faites pour faire joli, en analyse secondaire, mais que la conclusion ne dépend pas vraiment d’elles; on repose alors sur la présomption d’absence d’interaction (cf approche N°2) sécurisée par un choix pertinent de critères d’inclusion et d’exclusion (cf chapitre « les interactions qualitatives sont-elles fréquentes ? »). La non correction de multiplicité des tests est un moyen de remonter le risque alpha et d’améliorer la puissance, ce qui se rapproche du prior informatif (approche N°2). On pourrait aussi imaginer qu’on aille plus loin en remontant le seuil de significativité unilatéral de 2,5% à 5%, voire 10%. En pratique, le seuil de 5% bilatéral (2,5% unilatéral) est sacré; on ne peut pas y toucher; cependant on reste encore libre de tricher à mort avec la multiplicité des tests.

Approche N°2 : ajout d’une touche bayésienne

Cette approche N°2 n’est jamais employée même si elle pertinente; certainement en raison de l’obligation à intégrer une part de subjectivité. Cette seconde approche part du constat qu’avec de bons critères d’inclusion, les interactions qualitatives ou quantitatives fortes sont probablement rares. Pour illustrer cela, considérons le débat autour des vaccins anti-SARS-Cov-2 qui ont été tous principalement évalués chez des personnes jeunes (<= 65 ans) mais prescrits en premier chez des personnes âgés (>= 70 ans). D’abord, balayons une idée reçue ridicule. On ne peut pas conclure à l’efficacité d’un vaccin sur les personnes âgées en argumentant que 20% des patients étaient âgés et donc, qu’ils participent à l’analyse d’efficacité. En effet, si le vaccin était totalement inefficace chez les personnes âgées, alors les 80% de personnes jeunes pourraient très bien compenser cet effet nul du vaccin dans la statistique globale. Il est donc nécessaire de faire des analyses en sous-groupes, ou pas… De mon point de vue, subjectif, un vaccin efficace à plus de 90% chez des personnes d’âge moyen de 50 ans, pourrait être potentiellement un peu moins efficace sur des personnes âgés, mais une inefficacité totale chez la personne âgée paraît peu crédible; pas impossible, mais improbable. Je ne crois pas que ce phénomène d’inefficacité totale, passé un âge, ait jamais été rapporté pour un vaccin; les inefficacités majeures étant principalement rapportées chez des personnes fortement immunodéprimées. C’est ainsi, qu’avec une approche bayésienne, on aurait tendance à dire, qu’a priori, il y a de fortes chances que le vaccin ait une efficacité satisfaisante sur les personnes âgées, quand bien même l’échantillon n’en contient aucune ! L’expérience en population générale sur ces vaccins semble d’ailleurs avoir confirmé l’efficacité sur la personne âgée.

Une fois ce constat dressé, comment peut-on le traduire à termes statistiques ? Eh bien, on peut estimer les effets en sous-groupes avec une méthode bayésienne à prior informatif, considérant que la probabilité d’un effet proche de l’effet global est assez forte alors que la probabilité d’un effet nul ou négatif est faible, voire très faible. Le gros problème, c’est de régler la force de ce prior. Un prior d’informativité nulle, sera équivalent à l’approche N°1. En cas de résultat significatif, les conclusions seront très robustes, convainquant les bayésiens les plus sceptiques, mais bien souvent, les résultats ne seront pas concluants; ce qui n’est pas toujours un mal si on n’a réellement pas la donnée pour conclure. Au-delà une certaine force, le prior conduira à la conclusion d’un bénéfice en l’absence totale de donnée sur le sous-groupe; c’est-à-dire qu’on conclura au bénéfice du traitement avec N=0 observation. Ça peut paraître délirant, mais c’est loin de l’être; c’est ce qui a été fait par les autorités de santé publique qui ont recommandé les vaccins anti-SARS-Cov-2 aux personnes âgées alors que la quantité de données sur le sujet était négligeable ou nul. Un tel prior peut quand même conduire à remettre en cause l’efficacité du traitement si on dispose de suffisamment de données allant dans un sens contraire dans ce sous-groupe. En revenant au vaccin anti-SARS-Cov-2 et à la personne âgée : on part du principe que le vaccin est efficace chez la personne âgée, jusqu’à preuve du contraire; mais si on trouve une preuve du contraire, on est prêt à l’accepter. Même là, on peut régler le prior plus ou moins fortement, remettant en cause le bénéfice pour des preuves contraires faibles, ou au contraire exigeant des preuves contraires fortes.

Néanmoins je pense qu’il faut bien garder à l’esprit la distinction entre preuve et présomption. De mon point de vue, une conclusion à une efficacité basée sur N=0 est une présomption pas une preuve. Une présomption peut parfois être assez solide et des preuves peuvent être fragiles, mais la distinction pour moi, porte sur la source de données : la présomption est fondée sur un prior, et a une grand part de subjectivité alors que la preuve, même fragile, repose sur les données de l’étude (information contenue dans l’étude). Cette distinction est importante parce que la manière de remettre en cause la connaissance diffère. La preuve est remise en cause par la critique méthodologique de l’article alors que la présomption est remise en cause par la faible validité des connaissances externes ou surtout de leur articulation qui ont conduit à cette présomption.

C’est aussi sur cette présomption que tous les cliniciens reposent tous les jours en prescrivant des traitements à des patients pour lesquels il n’y a jamais eu d’analyse en sous-groupes ou pour lesquelles elles ne sont pas du tout concluantes ; c’est-à-dire la grande majorité des sous-groupes.

Comme cette approche N°2 n’est jamais appliquée formellement en pratique, que je ne l’ai moi-même jamais appliquée, je ne saurais pas vous conseiller sur la mise en oeuvre, et notamment le choix du prior.

Approche N°3 : échouer joyeusement au test d’interaction

Présentation

C’est cette approche N°3 que le reviewer m’a proposée. Elle viole allègrement les principes fréquentistes mais peut avoir une interprétation intéressante si on la regarde sous un autre angle. La procédure est la suivante:

Faire un test d’interaction, dont l’hypothèse nulle est : l’effet du traitement est identique dans tous les sous-groupes (p.e. identique chez les personnes âgées et les jeunes). Si on échoue à rejeter cette hypothèse nulle, alors l’accepter comme une vérité solide et conclure à l’équivalence des effets dans tous les sous-groupes. Cela permet alors de ne plus présenter les effets des sous-groupes mais juste un unique effet principal que l’on considère comme applicable à tous les sous-groupes. À l’opposé, si le test d’interaction est significatif, alors on considère que l’effet principal n’a plus de sens, et on présente seulement les effets des sous-groupes.

Problème conceptuel : acceptation de l’hypothèse nulle en situation de sous-puissance

Cette approche repose sur l’acceptation de l’hypothèse nulle dans une situation de puissance statistique catastrophique. On n’a pas les moyens de voir quoi que ce soit et on se rassure faussement en disant qu’on n’a rien vu. En fait, la puissance du test d’interaction est généralement très inférieure à celle des tests réalisés dans l’approche N°1.

Considérons par exemple une population dans laquelle il y a deux sous-groupes équilibrés (prévalences 50% / 50%), ce qui est un scenario très favorable. Une puissance à 90% aurait été calculée pour l’analyse principale, poolant les deux sous-groupes avec un effet +X. L’effet réel serait effectivement égal à +X dans le premier sous-groupe mais serait nul dans le second sous-groupe. L’approche N°1 conduirait à une puissance à 63% pour le premier sous-groupe:

pnorm((qnorm(0.975)+qnorm(0.90))/sqrt(2)-qnorm(0.975))

Avec l’approche N°1, le risque alpha serait à 5% pour le second sous-groupe et la notion de puissance ne s’appliquerait pas, puisque l’effet du sous-groupe est nul.

Avec l’approche N°3, la puissance du test d’interaction serait de 37%:

pnorm((qnorm(0.975)+qnorm(0.90))/sqrt(2*2)-qnorm(0.975))

En effet, ce test d’interaction est basé sur la différence d’effet entre deux groupes indépendants, avec les variances des effets qui se cumule; pour deux groupes de même taille avec des effets estimés de même précision, cela double la variance. La puissance diminue encore beaucoup si l’interaction est quantitative, avec, par exemple, un effet deux fois plus petit pour le second groupe que pour le premier. En effet, l’interaction à estimer vaudra alors X/2 plutôt que X, conduisant à une puissance à 12,5%:

pnorm((qnorm(0.975)+qnorm(0.90))/sqrt(2*2*4)-qnorm(0.975))

Et cela est pourtant estimé dans un scenario favorable, avec des groupes bien équilibrés, peu nombreux (2 sous-groupes seulement), et une puissance initialement calculée à 90%. Face à une puissance aussi ridicule, un test significatif perd en pertinence, puisque cette significativité est presque aussi probable sous l’hypothèse nulle que sous l’alternative; il n’est alors plus discriminant entre les deux et un résultat significatif ne fournit presque plus d’information sur la valeur de véracité de l’une ou l’autre des hypothèses. Il faut donc ignorer complètement les résultats d’un test significatif comme d’un test non significatif.

Interprétation abusive et comparatisme

Certes, l’interaction quantitative est beaucoup moins préoccupante que l’interaction qualitative, mais un test d’interaction échoué va donner l’illusion d’une preuve d’homogénéité des effets, comme le reviewer mentionnait, conduisant à généraliser l’intervalle de confiance de l’analyse poolée à tous les sous-groupes. Il me paraît plus honnête de se limiter à prouver qu’il existe un effet bénéfique dans chaque sous-groupe, sans les comparer les uns aux autres, ce qui est déjà très difficile à faire (approche N°1).

Ensuite, les tests d’interactions quantitatives ont d’autres problèmes. Sur de très grands échantillons, ce qui n’est pas courant mais arrive parfois, ils peuvent avoir une très bonne puissance et conclure à la non homogénéité de l’effet des traitements; rien de mal à ça, sauf que ça risque de conduire à la pratique du comparatisme; tendance à conclure, à tort, qu’un sous-groupe ne bénéficie pas du traitement seulement parce que son bénéfice est inférieur à celui d’un autre groupe. De mon point de vue, en situation de grand échantillon (aucune incertitude), l’évaluation de la balance bénéfices/risques dans un sous-groupe ne doit pas dépendre des autres sous-groupes. Si on considère qu’un bénéfice de +30% est bon, alors le sous-groupe A qui aura un effet de +31% aura une indication au traitement, peu importe que le groupe B ait un effet de +80%, de +31% ou de +15%. C’est l’approche N°1 qui permet de raisonner de cette manière alors que l’approche N°3 encourage très fortement à ne s’intéresser qu’à la comparaison des effets les uns aux autres.

Validité statistique conditionnelle

Dans cette approche, on n’évalue séparément les effets des différents groupes, seulement si le test d’interaction est significatif. En raison de la faible puissance qui y est généralement associée, le test ne sera significatif que si l’hétérogénéité observée est très largement supérieure à la l’hétérogénéité réelle (dans la population). L’interaction aura l’air énorme alors qu’il s’agira en réalité d’une interaction finalement modeste, voire nulle. Dans les cas où il existe une interaction réelle et modeste, mais la puissance est très faible (< 10%), il y aura de forts risques que l’interaction significative observée soit de signe opposé à l’interaction réelle et qu’on conclue, à tort que le sous-groupe de moindre bénéfice et celui de plus fort bénéfice !

Considérons une étude qui a été conçue avec une puissance à 80% pour mettre en évidence un effet +X. Cet effet s’avère être de +X dans le sous-groupe principal A représentant 75% de l’échantillon mais n’est que de 0.80×X dans un sous-groupe complémentaire B représentant 25% de l’échantillon. La probabilité d’observer une interaction significative du même signe que l’interaction réelle est de 4,3%:

pnorm((qnorm(0.975)+qnorm(0.80))/sqrt(1/0.25 + 1/0.75)*(1-0.80)-qnorm(0.975))

alors que la probabilité d’observer une interaction significative de signe opposé à l’interaction réelle est de 1,4%:

1-pnorm((qnorm(0.975)+qnorm(0.80))/sqrt(1/0.25 + 1/0.75)*(1-0.80)+qnorm(0.975))

Conditionnellement au fait que l’interaction soit significative, on a quand même 24% de risque de conclure à l’existence d’une interaction dans le sens opposé au sens réel de l’interaction !

Conditionnellement à un résultat significatif de signe correct (c’est-à-dire dans les 4,3% de bons scenarii), on estimera en moyenne l’effet du groupe A à 1.44×X (alors que c’est en réalité 1×X) et l’effet du groupe B à −0.51×X (alors que c’est en réalité 0.80×X) et on estimera la différence à 1.95×X en moyenne alors qu’elle est en réalité c’est 0.20×X. On conclura donc à une interaction majeure et on doutera très fortement du bénéfice dans le groupe B alors qu’en réalité il est presque équivalent à celui du groupe A. Il est même possible que certains concluent à tort à l’existence d’une interaction qualitative.

Conditionnellement à un résultat significatif de signe incorrect (dans 1,4% de cas), on estimera en moyenne l’effet du groupe A à 0.47×X (alors qu’en réalité c’est 1×X) et l’effet du groupe B à 2,38×X (alors qu’en réalité c’est 0.80×X), ce qui conduira à des estimations fortement erronnées mais heureusement ne conduira pas à la conclusion de l’existence d’une interaction qualitative en général.

L’histogramme ci-dessous est très parlant, car il montre que la directive du reviewer consistant à ne regarder l’interaction qu’en cas de significativité conduit à une estimation de ce terme d’interaction qui est tout sauf pertinente:

Alors que l’interaction réelle est de 0.20, l’interaction estimée sera toujours très grande ou très petite : à coup sûr inférieure à -1.62 ou supérieure à +1.62, mais jamais entre les deux.

Le code R qui permet de réaliser ces simulations est:

# simulations distributions des effets
# dans chacun des deux groupes
r1=rnorm(1e7, mean=1.00, sd=sqrt(1/0.75)/(qnorm(0.975)+qnorm(0.80)))
r2=rnorm(1e7, mean=0.80, sd=sqrt(1/0.25)/(qnorm(0.975)+qnorm(0.80)))
# standard error du terme d'interaction
SE=sqrt(1/0.75 + 1/0.25)/(qnorm(0.975)+qnorm(0.80))

ok=abs(r1-r2)/SE > qnorm(0.975) # significatif

# significatif dans le bon sens
round(mean(r1[ok & (r1-r2)>0]),2)
round(mean(r2[ok & (r1-r2)>0]),2)

# significatif dans le mauvais sens
round(mean(r1[ok & (r1-r2)<0]),2)
round(mean(r2[ok & (r1-r2)<0]),2)

# histogramme
hist(r1[ok]-r2[ok],breaks=200,prob=TRUE,
	las=1,main="",xlim=c(-4,4),
	xaxs="i",yaxs="i",ylim=c(0,2),
	ylab="Densité de proba",
	xlab="Distribution de l'interaction observée conditionnelle à sa significativité (gris)")
text(x=0.2,y=1.9,"Interaction réelle",adj=c(-0.1,0),col="blue")
abline(v=0.2,col="blue",xlab="Interaction observée",lwd=2)

J’ai aussi créé un code presque exact, basé sur la fonction de densité de proba de la loi normale, approchée par une loi discrète:

x=seq(-5, 5, 0.01)
g1=prop.table(dnorm(x, mean=1.00,
	sd=sqrt(1/0.75)/(qnorm(0.975)+qnorm(0.80))))
g2=prop.table(dnorm(x, mean=0.80,
	sd=sqrt(1/0.25)/(qnorm(0.975)+qnorm(0.80))))
SE=sqrt(1/0.75 + 1/0.25)/(qnorm(0.975)+qnorm(0.80))

df=cbind(expand.grid(x1=x,x2=x),
	expand.grid(p1=g1,p2=g2))
df$delta=df$x1-df$x2
df$signif=abs(df$delta)/SE >= qnorm(0.975)
df$proba=df$p1*df$p2
xdf=subset(df, signif)

ok=xdf$delta>0
round(sum(prop.table(xdf$proba[ok]) * xdf$x1[ok]),2)
round(sum(prop.table(xdf$proba[ok]) * xdf$x2[ok]),2)
round(sum(prop.table(xdf$proba[ok]) * xdf$delta[ok]),2)

ok=xdf$delta<0
round(sum(prop.table(xdf$proba[ok]) * xdf$x1[ok]),2)
round(sum(prop.table(xdf$proba[ok]) * xdf$x2[ok]),2)

max(xdf$delta[xdf$delta<0])
min(xdf$delta[xdf$delta>0])

Hypothèse nulle difficile à croire

L’hypothèse nulle d’homogénéité parfaite des effets est difficilement plausible pour les analyses en sous-groupes concernant les facteurs pronostics. En effet, selon que l’on utilise un modèle ou un autre, il y aura une interaction ou pas. Par exemple, si on considère un pourcentage de résultat favorable sous traitement à 30% dans le groupe A contre 20% sans traitement et un résultat favorable sous traitement à 40% dans le groupe B contre 30% sans traitement, les deux effets sont à +10% (modèle linéaire), mais les rapports de pourcentages, odds ratio ou effets dans un modèle probit, diffèrent. Ces modèles sont fondamentalement contradictoires sur leur conception de l’homogénéité des effets lorsque le niveau de base diffère selon le sous-groupe. Il n’y a aucune raison qu’un de ces modèles soit vrai car ils ont été simplement conçus pour des raisons de simplicité mathématique. En conséquence, l’absence totale d’interaction, dans l’un ou l’autre de ces modèles, paraît presque impossible. L’interaction peut être faible, mais elle ne devrait pas être nulle. C’est donc d’autant plus difficile d’accepter une hypothèse nulle d’homogénéité parfaite sachant que cette hypothèse est très peu plausible a priori.

Variante

Gail & Simon ont proposé un test d’interaction qualitative (Biometrics. 41 (2). 1985, pp. 361-372, doi:10.2307/2530862), qui évite la tendance au comparatisme et recentre le débat sur le problème des interactions vraiment préoccupantes. Malheureusement ce test aggrave énormément le problème de sous-puissance statistique et d’acceptation abusive de l’hypothèse nulle parce qu’il va globalement, ne fournir de résultat significatif que si un sous-groupe montre un effet significativement positif alors que l’autre montre un effet significativement négatif.

Synthèse des approches

L’approche bayésienne (approche N°2) est la plus générale et permet de bien comprendre les enjeux, même si elle est probablement difficile voire impossible à appliquer en pratique. L’approche N°1 est équivalente à une approche bayésienne (N°2) à prior non informatif. L’approche N°3 se rapproche de l’approche N°2 avec un prior absolument massif, complètement déraisonnable, et conduit à conclure abusivement à l’homogénéité de l’effet dans presque tous les cas; dans les rares cas où cette méthode arrive à montrer l’hétérogénéité de l’effet, les effets estimés sont extrêmement éloignés de la réalité, ce qui peut tromper encore plus que si on concluait à l’homogénéité. Dans la majorité des cas cette approche N°3 fait passer une présomption d’homogénéité pour une preuve, en renversant la charge de la preuve. Mais dans les rares cas où le test est significatif, la présomption (prior bayésien) s’efface brutalement conduisant à une conclusion complètement différente. C’est ainsi qu’une observation peut renverser brutalement (en passant de p=0.051 à p=0.049) les conclusions qui passent de « parfaite homogénéité » à « hétérogénéité absolument majeure », et cela, de manière très faiblement corrélée à l’existence d’une réelle hétérogénéité.

Notion d’effet moyen

Un problème avec l’approche N°3, c’est de ne pas assumer que même en présence d’interaction, l’effet poolé reste interprétable comme l’effet moyen du traitement donné à la population complète. Cet effet moyen est particulièrement pertinent lorsque le « traitement » est une exposition nocive ou bénéfique, difficile à contrôler individuellement, et qui, en conséquence est imposé à la population entière. À condition que l’échantillon soit représentatif, l’effet global représentera le bénéfice ou la nocivité collective. En cas de non représentativité de l’échantillon, on pourra éventuellement repondérer.

Pour bien interpréter un effet global comme un effet en sous-groupe, il est nécessaire de bien comprendre qu’en statistiques, on va toujours calculer des effets moyens, qui ne sont pas applicables à un niveau individuel. L’effet d’un traitement sera positif en moyenne s’il fait plus de bien que de mal, mais le bien et le mal ne concernent souvent pas les mêmes patients. Les effets indésirables sont très imprévisibles, de même que l’efficacité qui est très variable; c’est pour ça que même une étude en cross-over, effaçant les différences sur de très nombreuses variables observables comme inobservables entre les groupes, peut nécessiter de nombreux patients pour mettre en évidence une moyenne positive de l’effet. Philosophiquement, on doit pouvoir distinguer des effets « non déterministes ou apparentés » et les effets déterministes inconnus. Les effets « non déterministes ou apparentés » sont les effets que l’on peut considérer comme « totalement aléatoires », parce qu’ils ne seront jamais prévisibles, peu importe le nombre et la finesse des variables qu’on a mesurées. Après tout, le résultat d’un lancer de dé est totalement imprévisible, quand bien même la position et la vitesse initiale du dé sont parfaitement contrôlées et mesurées; même à un niveau subatomique, les mouvements d’un électron autour d’un noyau d’hydrogène sont totalement imprévisibles et doivent être modélisés par une fonction de densité de probabilité de présence (orbitale); le vivant est encore mois prévisible. Néanmoins, il existe aussi probablement des effets déterministes inconnus, sur des variables qui ne sont pas connues du tout ou pas mesurées habituellement (p.e. sous-type histologique de biologique moléculaire d’un cancer); il peut y avoir des interactions qualitatives sur ces variables, qui seront éventuellement connues à l’avenir. En l’absence de connaissance de ces variables d’interaction, on va évaluer l’effet du traitement sur un groupe comportant un mélange de patients avec un effet positif et de patients avec un effet négatif; selon le ratio des groupes et des effets, on pourra avoir un effet moyen positif. On considèrera alors qu’en moyenne, le traitement est bénéfique et on le donnera à tous les patients. Est-ce à tort ou à raison ? C’est une question philosophique à laquelle je répondrai qu’en l’absence de connaissance dur la variable inconnue, c’est à raison que l’on donne le traitement à tout le monde, car cela fournit la meilleure espérance au patient, compte-tenu des connaissances médicales actuelles; l’alternative est de ne pas donner le traitement, ce qui aura une espérance négative. La situation dans laquelle on a une puissance très insuffisante pour faire des analyses en sous-groupes n’est pas si différente; les connaissances actuelles ne permettent pas de répondre à la question du bénéfice dans chacun des sous-groupes, mais disent qu’en moyenne, il y a un bénéfice sur la population globale; on fait alors au mieux avec les connaissances actuelles. En réalité, ça dépend un peu de la force de la présomption d’interaction qualitative. Si cette présomption était forte, on s’obligerait à conclure sur les deux analyses en sous-groupes séparément (cf approche N°1). C’est pourquoi, par exemple, on ne mélange pas les enfants et les adultes dans l’évaluation des traitements psychiatriques, craignant fortement les interactions qualitatives.

Les interactions qualitatives sont-elles fréquentes ?

Les interactions qualitatives sont légion, mais on les voit rarement, parce qu’on les anticipe dans les critères d’inclusion et d’exclusion ! Depuis qu’on a abandonné la recherche de la panacée, on sait très bien que le traitement doit être adapté au malade parce qu’un traitement inadapté à de très forts risques de faire plus de mal que de bien. C’est pourquoi il est nécessaire de faire des diagnostics ! Selon l’étiologie d’une dyspnée, par exemple, le traitement sera très différent, et en administrant le mauvais traitement, on aura de forts risques d’aggraver l’état du patient. C’est pour ça qu’il faut dix ans d’étude pour avoir le droit de prescrire ces traitements : mal utilisés ils font plus de mal que de bien. Les projets de recherche prennent toujours en compte ces interactions en définissant, par des critères d’inclusion et de non-inclusion, une population suffisamment homogène pour que la présomption d’interaction soit faible. Dans les rares cas où la présomption d’interaction reste assez forte mais on espère néanmoins qu’il y aura un bénéfice dans deux sous-groupes, on peut planifier directement une analyse en sous-groupes selon l’approche N°1.

Quelques éléments supplémentaires sur le case report

Même si les exemples que j’ai pris sont calqués sur l’évaluation thérapeutique, la problématique est très semblable pour d’autres types d’étude, notamment les études diagnostiques. La question du reviewer portait justement sur une étude diagnostique dans laquelle on évalue l’apporte de l’intelligence artificielle au diagnostic radiologique de fractures. L’apport de l’intelligence artificielle, en sensibilité et spécificité, dépend potentiellement de la localisation de la fracture; le carpe ou les cotes ne sont pas le même problème, ni pour l’humain ni pour la machine. Contrairement à la plupart des études, celle-ci était très largement dimensionnée et conçue pour les analyses en sous-groupes. L’étude porte sur 480 examens radiographiques lus par 24 lecteurs. Tous les lecteurs lisaient tous les examens, avec et sans aide de l’intelligence artificielle avec une période de wash-out, ce qui conduisait à 24×480×2 = 23040 lectures. Comme la principale source de variance était le lecteur, les analyses en sous-groupes par localisation de l’examen radiographique portaient sur le même nombre de lecteurs (n=24) mais pas le même nombre d’examens (n < 480), ce qui baissait finalement modérément la puissance statistique par rapport à l’analyse principale. Au total, avec l’approche N°1, sur les 9 localisations considérées, on arrivait à montrer un apport positif de l’intelligence artificielle sur la sensibiltié pour 7 d’entre-elles et on échouait pour 2 d’entre elles au seuil de significativité 2,5% unilatéral (ou 5% bilatéral) sans correction de multiplicité des tests. En montant le niveau de significativité à 7,5% unilatéral (15% bilatéral), on montrait le bénéfice pour 8/9. Selon l’approche N°1, il est impossible de dire si pour les « 2 échecs », il n’y a pas de bénéfice de l’intelligence artificielle, s’il est plus faible que pour les autres localisations, ou si ce sont des échecs attribuables à de simples fluctutations d’échantillonnage défavorables. La présomption de forte interaction, voire d’une interaction qualitative, en raison de possibles effets pervers de l’aide, était loin d’être négligeable, mais on arrive quand même à avoir une conclusion assez forte sur 7 localisations et un doute pour les deux dernières.