Bon usage des nombres à virgule flottante : exemples

Calculer une variance peut paraître simple, mais il y a de très mauvaises formules qui circulent.

Considérons une variance portant sur des nombres énormes avec une très faible variance. Un cas fréquent, ce sont les dates ou paires dates/heures ! Ces valeurs sont typiquement représentées par un nombre de jours, secondes ou millisecondes depuis une certaine époque (1er janvier 1970, 1er janvier 1900 ou 1er janvier 1601). Le plus extrême est obtenu par les dates de création/modification des fichiers Windows qui sont enregistrés comme le nombre d’intervalles de 100 nanosecondes depuis le 1er janvier 1601. Ainsi, pour les dates récentes, les valeurs dépassent les 10 puissance 17.

Même si les dates sont très grandes, on peut vouloir les manipuler sur un intervalle de temps très court, tel que les quelques secondes d’une expérience de marche, comme j’ai pu voir dans un projet de recherche.

Considérons donc le code R suivant (exécuté sur R 4.0.3 sous Windows 10), basé sur deux dates/temps, exprimées en nombre de secondes depuis le 1er janvier 1970:

> x1=strptime("2021-01-01 00:00:01", format="%Y-%m-%d %H:%M:%S")
> x2=strptime("2021-01-01 00:00:02", format="%Y-%m-%d %H:%M:%S")
> var(as.numeric(c(x1,x2)))
[1] 0.5

Tout semble fonctionner bien avec la fonction var de base. Cependant si on réimplemente cette fonction avec la formule que l’on retrouve dans tous les livres, on obtient un résultat aberrant:

> badvar=function(x) (sum(x^2)-sum(x)^2/length(x))/(length(x)-1)
> badvar(as.numeric(c(x1,x2)))
[1] 0

Cela est explicable par le fait qu’on fait la soustraction de deux nombres très grand (sum(x^2) et sum(x)^2), dont la différence devrait être petite.

> sum(x^2) # représentation inexacte (dépasse 2^53, soit 9.0e15)
[1] 5.180695e+18
> sum(x)^2/(length(x))
[1] 5.180695e+18
> sum(x^2) - sum(x)^2/(length(x))
[1] 0

Une meilleure fonction commence par enlever la moyenne à toutes les valeurs, en utilisant la formule la plus naïve de la variance:

> bettervar=function(x) {sum((x-mean(x))^2)/(length(x)-1)}
> bettervar(x)
[1] 0.5

Sur de très grands vecteurs cette fonction peut quand même avoir une perte de précision parce que sa manière de calculer la moyenne des carrés, passe par une somme très grande.

Pour mieux comprendre ça, comparons l’implémentation de R de la fonction mean, à une implémentation naïve:

> badmean=function(v) {sum(v)/length(v)}
> v=rep(c(0,2^50 + 100), 2^24)
> mean(v[1:2])-2^49 # correct
[1] 50
> mean(v)-2^49 # correct aussi
[1] 50
> badmean(v)-2^49 # incorrect
[1] 7.5
> rowMeans(matrix(nrow=1,v))-2^49 # incorrect
[1] 7.5

En effet, la fonction sum() va conduire à un nombre aux alentours de 2^74, avec une représentation intermédiaire inexacte. On peut corriger ça, en rectifiant le biais, par deux appels à badmean.

> bettermean=function(v) {m = badmean(v); bias = badmean(v - m) ; m+bias}
> bettermean(v)-2^49
[1] 50

En effet, on commence par calculer une moyenne approximative avec badmean, puis on calcule la moyenne des écarts à cette moyenne, qui devrait correspondre à l’opposé du biais de la moyenne. Ce biais est estimé avec des sommes qui oscillent autour de zéro et ne vont donc pas dans des ordres de grandeur trop grands.

Sauf qu’il y a quand même un problème… Ces oscillations de la somme autour de zéro disparaissent si les valeurs sont plus ou moins ordonnées:

> w=sort(v)
> mean(w)-2^49 # incorrect
[1] 20.375
> bettermean(w)-2^49 # incorrect
[1] 20.375
> mean(w) - mean(v) # ouille
[1] -29.625

Ainsi, la fonction bettermean, comme la fonction mean de R, dépendent de l’ordre des valeurs dans le vecteur ! Heureusement sur plateforme PC, les FP80 utilisés par sum() fournissent une marge de sécurité durant le processus de sommation des nombres positifs et négatifs, de telle sorte que ce problème n’arrive que pour les grands vecteurs.

Par ailleurs la fonction badmean() dépend aussi de l’ordre des valeurs, comme dans l’exemple ci-dessous :

> v=rep(c(-2^52-1, 2^52), 1e6)
> badmean(v)
[1] -0.5
> badmean(sort(v))
[1] -0.002048

On peut alors implémenter encore une variance plus robuste en utilisant la fonction mean() deux fois:

> badvar=function(x) (sum(x^2) - sum(x)^2/length(x))/(length(x)-1)
> bettervar=function(x) {sum((x-mean(x))^2)/(length(x)-1)}
> bestvar = function(x) {mean((x-mean(x))^2)*(length(x)/(length(x)-1))}
> cst=2^25+10
> off=2^49
> v=rep(off+c(-cst, cst), each=1e7)
> vartheo = cst^2*(length(v)/(length(v)-1))
> bestvar(v) - vartheo # perfect
[1] 0
> bettervar(v) - vartheo # small error
[1] -87.5
> var(v) - vartheo # small error
[1] -87.5
> badvar(v) - vartheo # huge error
[1] -1.480078e+15
> badvar(v) # negative variance!
[1] -3.541775e+14

x87 Intel vs AMD

En développant divers tests de performances sur SSE2, AVX2 et les vieilles instructions x87, j’ai découvert un problème de performance avec les processeurs Intel. En présence de représentations spéciales (NaN, +Inf, -Inf), les processeurs Intel (Core i5-4460 de microarchitecture Haswell et Celeron J1900 de microarchitecture Atom/Bay Trail) verront un ralentissement de leurs calculs par un facteur 300 environ, alors que les processeurs AMD récents (Ryzen 1700X, Ryzen 3600) ne verront aucune dégradation de performances par rapport aux mêmes calculs sur des représentations non spéciales. Même si cela se voit très bien dans du code C, on peut aussi en voir l’impact avec du code R, notamment avec des NA.

> v=double(1e5) # sur Ryzen 1700 (Windows 10)
> system.time(for(i in 1:1e4)sum(v))
   user  system elapsed 
   1.51    0.00    1.51
> v[1]=NA;system.time(for(i in 1:1e4)sum(v))
   user  system elapsed 
    1.5     0.0     1.5 

Ainsi, il faut 1.5 secondes pour qu’un Ryzen 1700 fasse les 1e5*1e4 = 1 milliard d’additions ; peu importe la présence de NA. Cela ne montre pas le pic de performance FP80, en raison du vecteur de grande taille (800 Ko) qui dépasse de la cache L1 et L2 des processeurs modernes.

Le même code sur Haswell (Windows 10) montre des résultats très différents:

> v=double(1e5) # sur core i5-4460 (Haswell)
> system.time(for(i in 1:1e4)sum(v))
   user  system elapsed 
   0.93    0.00    0.95 
> v[1]=NA;system.time(for(i in 1:1e4)sum(v))
   user  system elapsed 
 134.04    0.17  134.33

Ainsi, les performances sont environ 144 fois plus basses dès qu’il y a un NA au début du vecteur ! Le même problème est observé sur Celeron J1900 (Ubuntu 64 bits).

> v=double(1e5) # sur Celeron J1900 (Atom/Bay Trail)
> system.time(for(i in 1:1e4)sum(v))
utilisateur     système      écoulé
      2.908       0.000       2.901
> v[1]=NA;system.time(for(i in 1:1e4)sum(v))
utilisateur     système      écoulé
    135.600       0.012     135.416
>

Ce problème de performance avec les NA est observé pour sum(), mean(), rowSums(), rowMeans(), prod(), mais pas pour var() et sd() qui doivent détecter les NA avant de lancer les calculs. Ces problèmes disparaissent dès qu’on utilise le paramètre na.rm=TRUE.

Je n’ai pas accès à un CPU Intel vraiment récent, mais je crains que ce problème soit toujours présent (si vous pouvez tester et répondre en commentaire, je serais très intéressé de savoir). On peut se demander si ça pourrait venir d’une mauvaise gestion des Exception Flags x87 par Intel.

Je ne suis pas le seul à avoir remarqué ça.

Piège des nombres à virgules flottante

J’ai l’habitude d’être très prudent avec la gestion des nombres à virgule flottante, notamment lorsqu’il y a des allers-retours avec des nombres entiers, mais je me suis récemment fait avoir comme un bleu sous R, probablement parce qu’une interruption de plusieurs jours dans mon écriture du script m’a fait oublier l’origine des données, ainsi qu’une autre problématique ayant détourné mon attention à un moment critique.

Afin de bien présenter cette problématique, je vais exposer tous les éléments en jeu en essayant de ne pas demander trop de pré-requis.

Représentation binaire d’un nombre entier

Commençons par décrire la représentation décimale des nombres entiers qu’on apprend à l’école primaire. Un nombre comme 532 comporte les trois chiffres « cinq », « trois » et « deux » dans le système décimal. Le nombre 532 est simplement égal à 5×100 + 3×10 + 2×1 = 5×10^2 + 3×10^1 + 2×10^0. De droite à gauche, les trois chiffres ont un « poids » croissant, selon les puissances de dix. Le chiffre 2 a un poids d’une unité, alors que le chiffre 3 représente les dizaines et le chiffre 5 représente les centaines. On peut représenter n’importe quel nombre entier positif ainsi. Le choix d’avoir dix chiffres (les chiffres de 0 à 9) plutôt que douze ou n’importe quel autre nombre est arbitraire. Une représentation bien plus simple pour les ordinateurs est la représentation binaire dans laquelle il n’y a que deux chiffres (0 et 1) et dans laquelle le poids des différents chiffres d’une représentation numérique correspondant à des puissances de deux. Ainsi, la représentation numérique binaire 11010 sera interprétée comme 1×2^4 + 1×2^3 + 0×2^2 + 1×2^1 + 0×2^0 = 1×16 + 1×8 + 0×4 + 1×2 + 0×0 = 16+8+2 = 26. Une autre représentation, assez populaire en informatique, est la représentation hexadécimale dans laquelle on a seize chiffres : 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F. Par exemple, la représentation numérique hexadécimale B16 correspondra au nombre entier 11×16^2 + 1×16^1 + 6×16^0 = 11×256 + 1×16 + 6×1 = 2838.

Pour la plupart de leurs calculs, les ordinateurs utiliseront, en interne, un format binaire avec un nombre fixé de chiffres. Par exemple, un ordinateur 32 bits manipulera majoritairement des nombres avec exactement 32 chiffres binaires, pas un de plus, ni un de moins. Le terme de bit est ici synonyme de chiffre binaire. Pour représenter les petits nombres, on utilisera des zéros comme préfixe. Par exemple, la valeur trois sera représentée par la séquence de 32 chiffres binaires suivante:

0000 0000 0000 0000 0000 0000 0000 0011

(les espaces ici sont seulement présentés pour faciliter la lecture mais ne sont pas signifiants)

En lisant cette représentation de droite à gauche, on peut calculer qu’elle représente le nombre 1×2^0 + 1×2^1 + 0×2^2 + 0×2^3 + … + 0×2^31 = 1 + 2 + 0 + 0 + … + 0 = 3

Le plus grand nombre que l’on puisse représenter ainsi ne contient que des 1 et est égal à la somme de toutes les puissances entières de deux, de 2^0 jusqu’à 2^31. On peut aisément montrer par récurrence que ce nombre est égal à 2^32 – 1, soit 4294967295 ; ce qui est suffisamment grand pour les décomptes de tous les jours. D’une manière générale, avec N bits, on peut représenter tous les nombre entiers de 0 jusqu’à 2^N-1, pour un total de 2^N nombres distincts. Il faut dire qu’une séquence ordonnée de N chiffres binaires dispose d’exactement 2^N combinaisons possibles (cela se démontre aisément par récurrence) si on considère que l’ordre a de l’importance, c’est-à-dire que 0011 n’est pas équivalent à 1100 ni à 1010. Pour les nombres entiers relatifs (positifs ou négatifs), il existe une représentation particulièrement astucieuse (représentation en complément à deux) si on connaît ℤ/nℤ, que je ne décrirai pas ici, car ce n’est pas nécessaire à mon exposé.

Un ordinateur moderne 64 bits, va principalement manipuler des nombres entiers 8 bits, 32 bits et 64 bits, selon le contexte.

Nombres à virgule flottante

Représentation

R, comme à peu près tous les logiciels et langages de programmation gèrent principalement les nombres réels de manière approximative sous forme de nombres à virgule flottante, communément appelés floating point values, ou FP. Cela permet de représenter approximativement des nombres comme 0.537877 ou 8.11215×10^56.

Les FP sont omniprésents en informatique : que ce soit sur une calculatrice, le processeur graphique (GPU) qui affiche un jeu vidéo sur mobile, PC, Mac ou console, n’importe quelle vidéo sur n’importe quelle plateforme postérieure à l’an 2000, ainsi que de nombreux formats de compression de sons ou d’image avec perte ; enfin, c’est encore utilisé par les réseaux neuronaux convolutionnels. C’est pourquoi, sur une vie entière au XXIème siècle, on est forcément confronté à un moment ou à un autre aux problèmes de ces représentations numériques. Enfin, un statisticien est amené à manipuler des FP dès qu’il utilise des variables quantitatives continues ou qu’il calcule des moyennes, pourcentages, odds ratio, hazard ratio ou toute autre statistique qui ne soit pas structurellement un nombre entier. Même si on peut utiliser des FP pendant des années sans en connaître les problèmes, on se cassera forcément une dent dessus à un moment où à un autre de sa carrière, parce que les FP sont structurellement imparfaits : ils tentent de résoudre au moins mal le problème insoluble de la représentation des nombres réels, appartenant à un ensemble infini indénombrable, dans un espace numérique fini et très borné.

Comme pour les nombres entiers, les représentations utiliseront typiquement des formats avec un nombre de bits bien prédéfini. R utilisera majoritairement des nombres à virgule flottante à double précision sur 64 bits, connus sous le nom de FP64 (conformes à la norme IEEE-754). Ces 64 bits sont divisés en trois partie : la mantisse, sur 52 bits, l’exposant sur 11 bits et le signe, sur 1 bit, pour un total de 52+11+1 = 64 bits. Avec ces trois composants, on représentera le nombre : signe × mantisse × 2^exposant. Le signe vaudra -1 ou +1. La mantisse sera un nombre entier positif. L’exposant sera un nombre entier positif ou négatif. Par exemple on pourra représenter -0,75 comme -1×3×2^-2. On constate, de prime abord, qu’il y a plusieurs représentations possibles pour le même nombre. On pourra représenter 16 avec une mantisse à 16 et un exposant à 0 (16×2^0), ou avec une mantisse à 8 et un exposant à 1 (16×2^1), ou encore, avec une mantisse à 4 et un exposant à 2 (4×2^2). Il existe un processus de normalisation, qui force l’usage d’une seule représentation avec une mantisse qui débute par le chiffre binaire 1 (sauf valeur trop petite, appelée dénormale). Ce bit à 1 est implicite dans la représentation binaire, ce qui permet d’avoir une précision de 53 bits sur la mantisse quand bien même elle n’est représentée que sur 52 bits. En bref, en dédoublonnant les représentations, on gagne en espace de représentation. C’est assez technique, mais ça ne change pas le fond du problème.

Problèmes (bugs insolubles) des nombres à virgules flottante

Cette représentation de nombre à virgule flottante souffre des mêmes forces et défaut que la représentation « scientifique » à quatre chiffres significatifs que vous avez peut-être apprise au lycée. Dans cette représentation, on va représenter n’importe quel nombre sous la forme A,BCD×10^EF où A, B, C, D, E et F représentent des chiffres de 0 à 9, sauf A qui vaudra au moins 1 ; sauf pour représenter zéro. Cette représentation est très précise, en absolu, pour les tous petits nombres. Par exemple, on pourra exprimer un volume de 1,567 microlitre comme 1,567×10^-6 litre ; ici, on a la précision du nanolitre. Cette même représentation sera beaucoup moins précise pour les gros nombres. Par exemple, mille deux-cent trois mètres cubes seront exprimés comme 1,203×10^6 litres, avec une précision du mètre cube. Si on veut ajouter quinze litres aux mille deux-cent trois mètres cubes, on va obtenir 1,203015×10^6 litres qui seront arrondis à 1,203×10^6 litres pour ne garder que les quatre chiffres significatifs selon le principe d’arrondi à la plus proche valeur représentable. Au final, le nombre est inchangé par l’addition des quinze litres, en raison de l’erreur d’arrondi finale. Ce problème d’arrondi touche aussi bien les nombres binaires à virgule flottante tels que les FP64. Même s’ils sont capables de représenter des nombres très proches de zéro, tels que 1/(2^1000), soit 2^-1000, et des nombres très grands, tels que 2^1000, l’addition de chiffres qui sont d’ordre de grandeur très différent, conduit à une perte de précision dommageable. Notamment, si on ajoute un nombre très grand à un nombre très petit, puis qu’on retranche le nombre très grand, alors le nombre très petit peu perdre énormément en précision, voire être annulé.

Ci-dessous, des exemples sous R:

> 1e9 + 42 - 1e9 # ça reste tolérable
[1] 42
> 1e17 + 42 - 1e17 # on commence à perdre
[1] 48
> 5e17 + 42 - 5e17 # ça fait mal
[1] 64
> 1e20 + 42 - 1e20 # on perd tout
[1] 0

L’ordre dans lequel on fait les calculs se met à avoir une importance !

> 5e20 - 5e20 + 42
[1] 42
> 1e20 + 42 - 1e20
[1] 0

En effet, la première ligne commence par le calcul (5e20 – 5e20) qui fait zéro, qui est ensuite additionné à 42, alors que la seconde ligne commence par 1e20+42 qui fait 1e20, puis retranche 1e20, ce qui conduit à la valeur zéro.

Un autre problème, c’est l’impossibilité de représenter de manière exacte certaines fractions pourtant simples, telles qu’un tiers. En effet, un tiers est approximé, en décimal, par 0.3333333333. Pour une représentation exacte, il faudrait une infinité de chiffres après la virgule. Même si 1/5 peut être représenté en décimal, il ne peut pas être représenté en binaire, parce que 5 est premier avec 2.

C’est ainsi, que des calculs très simples peuvent conduire à des résultats légèrement erronés:

> (1/5-1+1)*5 - 1
[1] -2.220446e-16

On peut se dire que le résultat est « suffisamment » proche de la bonne réponse (zéro) pour que ce soit tolérable, mais ça peut, en réalité, conduire à des réponses totalement erronées dès qu’on utilise l’opérateur d’égalité:

> (1/5-1+1)*5 == 1
[1] FALSE

Si R disposait d’un moteur de calcul exact, il afficherait TRUE ici. Sa réponse, d’une certaine manière, est complètement erronée.

Pour rendre les choses plus compliquées, R n’est pas honnête sur le vrai contenue des valeurs numériques, les arrondissant à quelques chiffres après la virgule (7 chiffres significatifs par défaut) pour l’affichage comme le montre le code ci-dessous:

> x=(1/5-1+1)*5
> x # on dirait qu'on a bien la valeur 1
[1] 1
> 1 == x # mais ce n'est pas vraiment 1
[1] FALSE
> x-1 # en fait, c'est un poil inférieur à 1
[1] -2.220446e-16
> 1-1 # alors que le "vrai" 1 est bien égal à 1
[1] 0

Ce comportement peut être corrigé par l’option digits, comme décrit ci-dessous:

> options(digits=22)
> (1/5-1+1)*5
[1] 0.99999999999999978
> 1/5
[1] 0.20000000000000001
> 1/5-1+1
[1] 0.19999999999999996

Là, on comprend beaucoup mieux les erreurs d’arrondi !

Il faut donc considérer que TOUT calcul avec nombre à virgule flottante est approximatif et que deux manières très légèrement différentes d’arriver au même calcul conduiront à deux approximations légèrement différentes. En fait, dans certains langages de programmation (p.e. le langage C), on n’a même pas la garantie qu’en réalisant deux fois exactement le même calcul on arrive au même résultat, parce que le compilateur peut traduire le code de manière différente.

Pour « tester » l’égalité, il faut toujours se donner des marges de manoeuvre, en comparant la valeur absolue de la différence à une valeur seuil correspondant à l’erreur maximale attendue.

Par exemple, plutôt que de tester (1/5-1+1)*5 == 1, on testera abs((1/5-1+1)*5 – 1)<1e-14.

On ne peut pas donner de valeur seuil universelle, parce que la bonne constante dépend de l’échelle sur laquelle on se trouve et de l’erreur cumulée des calculs. C’est-à-dire, que quand on manipule des nombres de l’ordre de 10^20, une erreur de plusieurs unités, voire plusieurs dizaines est rapidement arrivée, alors que quand on manipule correctement des nombres de l’ordre de 10^-20, alors on s’attend à avoir des erreurs inférieures à 10^-30. Parfois, on peut comparer le rapport entre les deux nombres à une valeur précise, mais ça ne marche que si les nombres restent écartés, d’une certaine manière, de la valeur zéro. Ce problème de connaissance de la bonne échelle se retrouve dans des problématiques apparentées, telles que la détection de convergence par stabilisation d’un paramètre (cf https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17885).

Il faut aussi savoir que même lorsqu’on manipule des nombres du même ordre de grandeur, les erreurs d’arrondi ont tendance à se cumuler. Même si les erreurs positives compensent en partie les erreurs négatives, grace à la méthode d’arrondi au plus proche, plus on fera d’opérations, plus l’erreur sur le résultat final sera importante.

Il existe une seule situation, fréquente dans R, pour laquelle les calculs peuvent être considérés comme exacts : lorsqu’on manipule des nombres entiers suffisamment petits (< 2^52) sous forme de nombre à virgule flottante, en veillant à ne faire aucune opération qui les fasse passer en nombre non entiers. En effet, dans R, si on oublie de préciser un suffixe ‘L’ derrière une constante numérique, elle représente un nombre à virgule flottante.

> storage.mode(42L) # vrai nombre entier
[1] "integer"
> storage.mode(42) # nombre à virgule flottante
[1] "double"

La manipulation du nombre à virgule flottante 42 est relativement sûre. J’appelerai ces nombres à virgule flottante représentant des entiers, des pseudo-entiers.

> 42*2 - 42 == 42 # ça reste sûr
[1] TRUE

Les opérations acceptables avec les pseudo-entiers sont : l’addition, la soustraction, la multiplication. L’exponentiation par l’opérateur ^ est acceptable si on reste dans des valeurs entières pas trop grandes (en restant en dessous de 2^52 pour le résultat final) et qu’on reste sur plateforme PC. En effet, elle repose sur la fonction pow() en langage C dont je crains qu’elle soit fortement approximative sur beaucoup de plateformes, en passant par des fonctions logarithme et exponentielles qui n’aient pas toute la précision des FP64. Cet opérateur ^ est toujours risqué car il n’existe pas en version vraiment entière mais seulement sur les nombres à virgule flottante, comme le montre le code suivant:

> storage.mode(3L^3L)
[1] "double"

Les pseudo-entiers ne marchent plus correctement à partir de 2^54:

> x=18014398509481984
> x+1 == x
[1] TRUE

R fournit quelques points de repères sur la précision des nombres à virgule flottantes dans l’objet .Machine.

Par exemple la constante .Machine$double.eps environ égale à 2.22e-16 est la plus petite valeur x telle que 1 + x != 1 selon la documentation. En réalité, ce n’est pas tout à fait vrai, il s’agirait plutôt de la plus petite valeur telle que 1 + x -1 == x.

Représentations spéciales

Outre les problèmes d’arrondi, les FP ont aussi des problèmes d’overflow et underflow parce que l’exposant est limité (variant de -1023 à +1024). La valeur 1,8×10^308 est le plus grand nombre représentable sur FP64 alors que 5×10^-324 est le nombre le plus petit (proche de zéro) représentable. Si on dépasse le plus grand nombre, on fait un « overflow » qui conduit à une valeur spéciale Inf représentant l’infini. Si on passe en dessous du plus petit nombre représentable, alors on atteint la valeur zéro par un « underflow ».

> 1e308*100
[1] Inf
> 1e-324/100 == 0
[1] TRUE

Il existe deux zéros différents et deux infinis différents : positifs et négatifs, avec des relations assez attendues entre les quatre:

> 1/0
[1] Inf
> 1/-0
[1] -Inf
> 1/Inf
[1] 0
> 1/(1/Inf)
[1] Inf
> 1/-Inf
[1] 0
> 1/(1/-Inf)
[1] -Inf

Néanmoins, les deux zéros sont égaux sans que leurs inverses le soient !

> pos=0
> neg=-0
> pos == neg
[1] TRUE
> 1/pos == 1/neg
[1] FALSE
> sign(pos)
[1] 0
> sign(neg)
[1] 0

La valeur Inf représente donc un nombre supérieur à tout ce qui est représentable mais inconnu. Alors que le +0 représente un nombre positif inférieur à tout ce qui est représentable mais inconnu. Certains calculs, tels que Inf – Inf conduisent à un résultat totalement inconnu qui est représenté par une valeur spéciale NaN. Cette valeur NaN se lit « Not a Number ».

> Inf - Inf
[1] NaN

Ce NaN se propage un peu de la même manière que les données manquantes (NA). La fonction is.na reconnaît d’ailleurs ces NaN comme des données manquantes. La fonction is.nan, par contre est spécifique des NaN.

Cas rapporté

Revenons à mon cas, qui illustre bien la perversité des FP si on n’est pas attentif. J’avais plusieurs fichiers de données correspondant à la capture de mouvement de plusieurs mires sur un mouvement de marche de plusieurs sujets. La capture de mouvement était faite cent fois par seconde par des caméras. Le timing de certains événements était enregistré en secondes, avec deux chiffres après la virgule. Ainsi, le temps 1,75 seconde correspondait à la frame numéro 175 de l’enregistrement. Une simple multiplication par 100 permettait de passer du timing enregistré au numéro de frame.

Pour un des calculs, j’avais besoin de calculer, sur un intervalle temporelle, la différence entre la coordonnée Z (hauteur) minimale observée sur l’intervalle [T1 ; T2] et la coordonnée Z observée à la fin de l’intervalle, c’est-à-dire au temps T2. Étant donné que T2 fait partie de l’intervalle, cette différence était obligatoirement positive ou nulle. J’observais pourtant des valeurs négatives ! Cela me conduit à un bout de code ressemblant à peu près à :

> min(Z[T1:T2])
[1] 23.78
> Z[T2]
[1] 21.14

(j’ai changé les noms des variables et les valeurs numériques, mais l’idée reste là)

Voilà à peu près ce qu’ont donné les investigations:

> T1 # ok
[1] 10
> T2 # ok
[1] 15
> T1:T2 # semble ok
[1] 10 11 12 13 14 15
> Z=1:15 # ça va simplifier l'interprétation
> Z[T2] # toujours ok
[1] 15
> Z[T1:T2] # tout est décalé vers le bas !
[1]  9 10 11 12 13 14
> Z[10:15]# mais là ça marche 
[1] 10 11 12 13 14 15

Pouvez-vous deviner le problème ?

Il vient de T1 et T2 qui ne sont pas vraiment entiers mais flottant, très proches de nombres entiers, d’une manière particulièrement perverse:

> T1-10
[1] -1.065814e-14
> T2-15
[1] 1.065814e-14

Et voilà !

Avec options(digits=22) on comprend encore mieux ce qui se passe:

> options(digits=22)
> T1
[1] 9.9999999999999893
> T2
[1] 15.000000000000011
> T1:T2
[1]  9.9999999999999893 10.9999999999999893 11.9999999999999893
[4] 12.9999999999999893 13.9999999999999893 14.9999999999999893
> as.integer(T2)
[1] 15
> as.integer(T1:T2)
[1]  9 10 11 12 13 14

Il faut comprendre que l’opérateur deux-points est un raccourci pour la fonction seq(from,to,by) qui marche avec des nombres non entiers. C’est ainsi, qu’on peut générer des séquences de nombres, en additionnant répététivement by à from, sans dépasser to. Ainsi, seq(0, 1, 0.40) va générer la sequence c(0, 0.40, 0.80). Eh bien, A:B est équivalent à seq(A, B, 1). Ici, le décalage vers le bas de T1 est répercuté sur la séquence entière. Ensuite, lorsqu’on indexe un vecteur par des indices non entiers, ils sont implicitement convertis par la fonction as.integer() qui tronque la virgule, c’est-à-dire, qui arrondit vers le bas pour les nombres positifs et vers le haut pour les nombres négatifs.

Le problème provenait des timings T1 et T2 que j’avais oublié d’arrondir au plus proche après la multiplication par 100 pour passer des secondes aux centisecondes.

Une fois détecté, le problème est trivial à résoudre, avec un simple appel à la fonction round().

Conséquences

Cette erreur a eu des conséquences non négligeables. J’ai rendu un rapport d’analyse statistique préliminaire dont les résultats étaient très légèrement erronés, car ils étaient fondés sur des valeurs qui avaient été calculés avec des décalages accidentels d’un centième de seconde. Ce n’est que lorsque j’ai voulu calculer des statistiques post hoc, que j’ai découvert l’incohérence ; celle-ci m’a ensuite conduit à découvrir le bug. Heureusement, une fois l’incohérence découverte l’investigation du bug n’a demandé que quelques minutes.

Ce qu’il faut en retenir

Il faut bien connaître les FP, parce qu’on va forcément se les prendre dans la face un jour ou l’autre. Je conseille aussi d’avoir une approche très dichotomique et systématique au débogage. Non seulement, cela permet de corriger son bug, mais parfois aussi d’en découvrir dans le logiciel et de faire un beau bug report !

Est-ce un bug de R ? Certains diraient que it works as designed ce à que je répondrais it is badly designed. Plutôt qu’un bug report, cela mériterait une feature request, pour demander au beau message d’erreur pour des indiçages aussi foireux.

Pour aller plus loin

Je conseille le site Web suivant: https://randomascii.wordpress.com/2012/04/21/exceptional-floating-point/

Update (14/05/2021)

Le logiciel R utilise pour certaines routines, sur plateforme PC (i386 ou x86_64), les FP80 de précision étendue offerts par l’interface x87 des microprocesseurs des PC plutôt que l’interface SSE2 qui implémente des FP64. C’est ainsi, que certaines routines fonctionnent avec une précision excédentaire.

> x=0;for(i in 1:1e9) {x=x+1/3};print(x); # FP64
[1] 333333332.6651181
> sum(rep(1/3,1e9)) # FP80 en interne
[1] 333333333.33365959
> 1/3*1e9 # pas de cumul d'erreur (en relatif)
[1] 333333333.33333331

Avec les FP80, l’erreur cumulée est 2048 fois plus faible, révélant une mantisse 63 bits plutôt que 52 bits, soit 11 bits supplémentaires, conduisant à une erreur 2^11 = 2048 fois plus faible.

Les fonctions rowSums, colSums, rowMeans, colMeans fonctionnent aussi en FP80 sur plateforme PC. On peut imaginer qu’on gagnerait nettement en performance si R utilisait des FP64 avec du SSE2 (SIMD 2×FP64) ou de l’AVX2 (SIMD 4×FP64), voire de l’AVX-512 (SIMD 8×FP64), mais en réalité, même pour un logiciel numérique comme R, le gain n’en vaut pas la chandelle par rapport au gain de précision fourni par les vieux FP80 des années 80 qui ont fait la renommée du Pentium. Pour ceux qui voudraient absolument booster les performances au détriment de la précision, il existe rowsums dans le package Rfast.

> Rfast::rowsums(matrix(nrow=1,rep(1/3,1e9)))
[1] 333333332.6651181

Mais le gain de performances est loin d’être extraordinaire sur des énormes vecteurs qui débordent largement de la cache L3, même lorsque la préparation de la matrice est minimaliste (testé sur AMD Ryzen 1700 qui dispose de 2 unités 128 bits capables de faire des additions 2×FP64):

> system.time(Rfast::colsums(matrix(nrow=1e9,ncol=1,1/3)))
   user  system elapsed 
   1.58    0.41    1.99 
> system.time(colSums(matrix(nrow=1e9,ncol=1,1/3)))
   user  system elapsed 
   2.56    0.49    3.05

Sur des vecteurs de taille modeste, tenant en cache L1, on peut perdre en performance:

> system.time(replicate(1e5,Rfast::colsums(matrix(nrow=1e3,ncol=1,1/3))))
   user  system elapsed 
   1.53    0.00    1.54 
> system.time(replicate(1e5,colSums(matrix(nrow=1e3,ncol=1,1/3))))
   user  system elapsed 
   1.13    0.00    1.13 

C’est probablement dû à l’inertie de lancement de la fonction Rfast::colsums. Cela nous montre que même sur des microbenchmarks, le gain SIMD est loin d’être extraordinaire.

Si on veut gagner un facteur deux en performances, on est obligé de créer un scenario spécifiquement optimisé pour ça:

> system.time(replicate(1e4,Rfast::colsums(matrix(nrow=1e5,ncol=1,1/3))))
   user  system elapsed 
   1.26    0.00    1.27 
> system.time(replicate(1e4,colSums(matrix(nrow=1e5,ncol=1,1/3))))
   user  system elapsed 
   2.33    0.00    2.33

Sur microarchitecture Intel, le gain devrait être nettement supérieur (https://www.agner.org/optimize/blog/read.php?i=838) depuis Haswell grace aux deux units 256 bits capables toutes deux de faire des additions de 4 FP64.

Update 2 (14/05/2021)

La fonction Rfast::colsums semble ne pas utiliser AVX2. Afin de vraiment évaluer les performances d’AVX2, j’ai dû développer un petit package:

#include <Rinternals.h>

#define NBYTES 32
#define NFP64 (NBYTES/8)
typedef double multifp __attribute__((vector_size (NBYTES)));

static double fpaddpar(double  *v, R_xlen_t n) {
        multifp psum1={0};
        multifp psum2={0};
        R_xlen_t i=0;
        double asum=0;
        
        /* ensure alignment */
        size_t uv = (size_t)v;
        while ((uv % NBYTES)!=0) {
                asum += *v++;
                uv = (size_t)v;
                n--;
        }

        multifp * restrict w=(multifp * restrict)v;
        R_xlen_t n2 = n/NFP64;
        for (i=0; i <= n2-2; i+=2) {
                psum1 += w[i];
                psum2 += w[i+1];
        }
                
        psum1 += psum2;
        #if NBYTES == 16
        	asum += psum1[0] + psum1[1];
        #elif NBYTES == 32
        	asum += psum1[0] + psum1[1] + psum1[2] + psum1[3];
        #else
        	#error number of bytes unsupported
        #endif
        
        for(i=i*NFP64; i < n; i++) {
                asum += v[i];
        }

        return asum;
}

static double slow_fpaddpar(double  *v, R_xlen_t n) {
  long double out=0;
  for(R_xlen_t i=0;i<n;i++) {
   out += v[i];
  }
  return out;
}

static SEXP fastsum(SEXP v) {
  SEXP result = PROTECT(allocVector(REALSXP, 1));
  REAL(result)[0] = fpaddpar(REAL(v), XLENGTH(v));

  UNPROTECT(1);
  return result;
}
static SEXP slowsum(SEXP v) {
  SEXP result = PROTECT(allocVector(REALSXP, 1));
  REAL(result)[0] = slow_fpaddpar(REAL(v), XLENGTH(v));

  UNPROTECT(1);
  return result;
}
static const R_CallMethodDef callMethods[]  = {
  {"fastsum", (DL_FUNC) &fastsum, 1},
  {"slowsum", (DL_FUNC) &slowsum, 1},
  {NULL, NULL, 0}
};

__declspec(dllexport) void R_init_fastsum(DllInfo *info) {
  R_registerRoutines(info, NULL, callMethods, NULL, NULL);
}

Le coeur de l’algorithme peut être retrouvé avec objdump -S pour avoir le code assembleur (usage de GCC 8.3.0 64 bits avec les options -O2 -march=haswell -mtune=haswell) :

 120:   48 83 c2 02             add    $0x2,%rdx
 124:   c4 c1 7d 58 00          vaddpd (%r8),%ymm0,%ymm0
 129:   c4 c1 75 58 48 20       vaddpd 0x20(%r8),%ymm1,%ymm1
 12f:   49 83 c0 40             add    $0x40,%r8
 133:   48 39 d1                cmp    %rdx,%rcx
 136:   7f e8                   jg     120 <fastsum+0x90>

Ainsi, AVX2 est bien utilisé, avec une parallélisation qui permet d’exploiter à fond les deux unités 256 bits du Haswell.

En créant des conditions artificielles absolument idéales, on voit effectivement un gain:

> v=rep(1/3,1e5) # AMD Ryzen 1700 à 3.00 Ghz
> system.time(replicate(1e4,.Call(C_fastsum,v)))
   user  system elapsed 
   0.14    0.00    0.14 
> system.time(replicate(1e4,sum(v)))
   user  system elapsed 
   1.61    0.00    1.61 

Mais dès qu’on prend en compte le temps de génération du vecteur, même simple, ce bénéfice s’estompe énormément:

> system.time(replicate(1e4,.Call(C_fastsum,rep(1/3, 1e5))))
   user  system elapsed 
   2.53    0.94    3.47 
> system.time(replicate(1e4,sum(rep(1/3, 1e5))))
   user  system elapsed 
   3.67    0.89    4.56

De même sur Haswell (core i5-4460 à 3.20 Ghz)

> v=rep(1/3,1e5)
> system.time(replicate(1e4,.Call(C_fastsum,v)))
   user  system elapsed 
   0.21    0.00    0.21 
> system.time(replicate(1e4,sum(v)))
   user  system elapsed 
   0.91    0.00    0.91 
> system.time(replicate(1e4,.Call(C_fastsum,rep(1/3,1e5))))
   user  system elapsed 
   2.28    1.30    3.58 
> system.time(replicate(1e4,sum(rep(1/3,1e5))))
   user  system elapsed 
   2.93    1.26    4.22

Étonnamment, le vieux Haswell core i5-4460 de 2014 est plus performant en FP80 que le Ryzen de 2017. En tout cas, les gains du SIMD se sentent toujours sur les microbenchmarks mais sont peu susceptibles d’avoir une influence dans du code R réel.

Update 3 (14/05/2021)

Après un déroulage de la boucle des FP80, la différence de performance s’amenuise:


#define NBYTES 32
#define NFP64 (NBYTES/8)
typedef double multifp __attribute__((vector_size (NBYTES)));

static double fpaddpar(double  *v, R_xlen_t n) {
        multifp psum1={0};
        multifp psum2={0};
        multifp psum3={0};
        multifp psum4={0};
        R_xlen_t i=0;
        double asum=0;
        
        /* ensure alignment */
        size_t uv = (size_t)v;
        while ((uv % NBYTES)!=0) {
                asum += *v++;
                uv = (size_t)v;
                n--;
        }

        multifp * restrict w=(multifp * restrict)v;
        R_xlen_t n2 = n/NFP64;
        for (i=0; i <= n2-4; i+=4) {
                psum1 += w[i];
                psum2 += w[i+1];
                psum3 += w[i+2];
                psum4 += w[i+3];
        }
                
        psum1 += psum2 + psum3 + psum4;
        #if NBYTES == 16
        	asum += psum1[0] + psum1[1];
        #elif NBYTES == 32
        	asum += psum1[0] + psum1[1] + psum1[2] + psum1[3];
        #else
        	#error number of bytes unsupported
        #endif
        
        for(i=i*NFP64; i < n; i++) {
                asum += v[i];
        }

        return asum;
}

static double slow_fpaddpar(double  *v, R_xlen_t n) {
  long double out1=0;
  long double out2=0;
  long double out3=0;
  long double out4=0;
  long double out5=0;
  long double out6=0;

  R_xlen_t i=0;
  for(;i<n;i+=6) {
   out1 += v[i];
   out2 += v[i+1];
   out3 += v[i+2];
   out4 += v[i+3];
   out5 += v[i+4];
   out6 += v[i+5];
  }
  out1+=out2+out3+out4+out5+out6;
  for(;i<n;i++) {
    out1+=v[i];
  }
  return out1;
}

Sur le core i5-4460 à 3.2 Ghz:

> v=rep(1/3, 1e5)
> system.time(replicate(1e4,.Call(C_fastsum,v)))
   user  system elapsed 
   0.22    0.00    0.22 
> system.time(replicate(1e4,.Call(C_slowsum,v)))
   user  system elapsed 
   0.42    0.02    0.44

Et sur le Ryzen 1700 à 3.0 Ghz (qui monte en réalité à 3.2 Ghz bien ventilé) :

> v=rep(1/3, 1e5)
> system.time(replicate(1e4,.Call(C_fastsum,v)))
   user  system elapsed 
   0.14    0.00    0.14 
> system.time(replicate(1e4,.Call(C_slowsum,v)))
   user  system elapsed 
   0.33    0.00    0.33 

Après analyse plus fine, les unités SIMD 128 bits (Ryzen) et 256 bits (Haswell) n’arrivent pas à saturation à cause du temps de load/store dans la cache L1/L2. Si on veut artificiellement montrer la supériorité du SIMD, il faut bidouiller le micro-benchmark pour que tout tienne de justesse en cache L1 et minimiser strictement les overheads d’appels en s’aidant de l’ALTREP de R 4.0 avec la boucle for.

> v=rep(1/3, 3e3) # Ryzen 1700
> A_fastsum=C_fastsum$address
> A_slowsum = C_slowsum$address
> system.time(for(i in 1:1e6) .Call(A_fastsum,v))
   user  system elapsed 
   0.33    0.00    0.33 
> system.time(for(i in 1:1e6) .Call(A_slowsum,v))
   user  system elapsed 
   1.08    0.00    1.08 

Même en trichant au maximum, on voit que c’est difficile de montrer un gros bénéfice du SIMD. L’architecture Zen étant limitée à 2 loads 128 bits/cycle depuis la cache L1 vers un registre.

Conditions de validité des estimateurs

J’avais vivement critiqué l’usage des tests de normalité dans un billet précédent, mais je n’avais pas insisté sur les bonnes manières de faire du cas de mésusage le plus fréquent des tests de normalité : s’assurer de la validité des estimateurs statistiques reposant sur la normalité asymptotique, tel que le test de Student sur séries appariées. Le test de normalité, dans ce contexte, fait exactement le contraire de ce qu’il faut, puisqu’il a une chance d’autant plus grande de conclure à la non normalité que l’échantillon est grand alors que la moyenne s’approche d’autant plus d’une loi normale. Ici, ce qui nous intéresse, c’est le biais d’estimation, exprimable par :

  1. Les risques de sous-estimation et surestimation de la statistique d’intérêt par l’intervalle de confiance
  2. Le risque alpha réel, qui doit s’approcher autant que possible du risque alpha nominal (généralement 5%)

Enjeu du choix d’un estimateur

Un point important est de bien distinguer l’estimateur du paramètre de la population qu’il estime. L’estimateur est une procédure qui, à partir d’un échantillon, fournit une estimation du paramètre, plus ou moins proche de cette cible. Dans les estimateurs, on distingue aussi:

  1. Les estimateurs ponctuels, qui fournissent une seule valeur, que l’on espère être aussi proche que possible du paramètre
  2. Les estimateurs d’intervalles de confiance, qui fournissent un intervalle qui doit avoir une chance élevée (p.e. 95%) de contenir le paramètre

Dans un projet de recherche on voudra généralement estimer un paramètre précis, mais le statisticien aura la liberté de choisir les estimateurs (p.e. ponctuel et intervalle de confiance) qui lui conviennent dans le paragraphe statistique du protocole de recherche. Par exemple, on voudra estimer la différence de moyennes entre deux groupes, mais le statisticien pourra calculer l’intervalle de confiance de cette différence de moyennes par:

  1. Test de Student sur séries indépendantes, avec hypothèse de variances égales
  2. Test d’Aspin-Welch qui est la variante à « variances inégales » du test de Student sur séries indépendantes
  3. Bootstrap non paramétrique, avec ses variantes : percentile, BCa, avec approximation normale, etc.
  4. Encore d’autres méthodes sont possibles, notamment dans des situations où on a des informations supplémentaires sur les distributions de ces moyennes, telles que le fait qu’elles suivent des lois binomiales ou de Poisson

Tous ces estimateurs estiment la même chose : différence de deux moyennes. Ils répondent à la même question de recherche et leur choix doit être guidé principalement par deux choses : les biais et l’efficacité. L’estimateur d’intervalle de confiance à 95% idéal conduira à une faible largeur d’intervalle de confiance tout en contrôlant parfaitement les risques de sous-estimation et surestimation à 2,5% chacun.

Choisir un estimateur ou un paramètre ?

Il faut bien distinguer le choix d’estimateur du choix du paramètre de la population à estimer. Le choix du paramètre de la population à estimer doit se faire en premier et il doit être guidé par la question de recherche. Le choix de l’estimateur doit se faire sur les biais et l’efficacité. Dans certaines situtations, le paramètre d’intérêt n’est tout simplement pas estimable car l’échantillon est trop petit (ou pour d’autres raisons). Au pire, il faut complètement remettre en cause la faisabilité de l’étude ; changer complètement de méthodologie, voire renoncer complètement à répondre à la question de recherche plutôt que de répondre de manière correcte à une question hors sujet.

Je vais donner un exemple caricatural qui permet de comprendre l’enjeu. Considérons une assurance souhaitant calculer la rentabilité d’une assurance portant sur des satellites artificiels face au risque de défaillance du lanceur. La rentabilité de l’assurance repose sur le fait que les primes, en moyenne, sont plus grandes que le coût des sinistres assurés plus les frais de fonctionnement. L’incertitude principale concerne la fréquence des sinistres. Si l’échantillon ne contient qu’une centaine d’événements assurés et aucun sinistre, il est impossible d’évaluer la rentabilité de l’assurance de manière empirique directe. Le paramètre à estimer est une différence de moyennes. Certains pourraient être tentés d’estimer une rentabilité médiane plutôt que moyenne, mais celle-ci serait hors sujet. Il est certain que moins de 50% des assurés auront un sinistre. Une prime d’assurance égale aux frais de gestion du dossier plus trois dollars aura l’air alors rentable.

Il va falloir parfois estimer la fréquence et le coût des sinistres de manière indirecte, en se basant sur des statistiques générales de taux de défaillance des lanceurs.

Validité a posteriori ou a priori

Pour commencer, il faut distinguer deux approches distinctes à la validité : approche a priori et approche a posteriori. L’approche bien trop souvent choisie est purement a posteriori dans laquelle les données observées sur un échantillon servent à calculer une ou plusieurs statistiques de validité. Si toutes les statistiques de validité sont satisfaisantes, selon des seuils prédéfinis ou selon une évaluation subjective, alors, l’estimateur statistique (p.e. Student) est considéré comme valide. Sinon, l’estimateur est considéré comme invalide et un autre estimateur statistique est utilisé, voire dans le pire des cas, le statisticien renonce à fournir une estimation.

Cette approche a posteriori est particulièrement problématique par rapport à la théorie fréquentiste qui doit juger de la procédure statistique dans son ensemble. Considérons par exemple, la procédure statistique suivante:

  1. Réaliser un test d’égalité des variances entre deux groupes
  2. Si les variances sont significativement différentes, alors réaliser estimer la différence de moyennes par Aspin-Welch
  3. Sinon, estimer la différence de moyennes par Student

Cette procédure statistique devient un nouvel estimateur dont on doit étudier les propriétés statistiques dans l’ensemble. Si vous étudiez cet estimateur sur des simulations, vous devriez constater qu’il est fortement biaisé en cas d’hétéroscédasticité forte sur de petits échantillons de taille déséquilibrée, parce que le test d’égalité des variances sera de puissance insuffisante, conduisant à une procédure proche d’un Student systématique.

Renoncer à fournir une statistique en raison d’une invalidité statistique fait aussi partie de la procédure dans son ensemble et engendre un selective reporting biais. Je pense notamment aux odds ratio dans une étude cas-témoin dont l’estimateur semble d’autant moins valide que l’odds ratio observé est éloigné de un. Ainsi, il existera une covariance entre l’odds ratio observé et ses chances de publication. Par exemple, considérons une étude cas-témoin avec 100 cas, 100 témoins et 20 exposés. Si on observe 10 exposés chez les cas et 10 exposés chez les témoins (odds ratio = 1), beaucoup de procédures considérerons que l’estimation de l’odds ratio est valide. Si on observe les 19 exposés chez les cas et 1 exposé chez les témoins, on peut craindre une mauvaise fiabilité de l’estimateur de Wald de l’intervalle de confiance dans une régression logistique. Si on rennonce à fournir cet odds ratio, cela est alors entièrement dû à la valeur observée trop forte de cet odds ratio ; en faisant ainsi, on biaise la méta-analyse, vers une sous-estimation de l’odds ratio réel.

On peut limiter ce problème en particulier, en basant le seuil décisionnel sur le nombre total de cas exposés (ici 20 exposés) plutôt que sur le nombre d’exposés dans les groupes de cas (n=19) et de témoins (n=1) ; cela évite que le selective reporting ait une covariance avec l’odds ratio observé et engendre un biais important dans la méta-analyse. Néanmoins, ça ne résoud pas complètement le problème de validité si jamais on obtient 19 exposés dans un groupe et 1 dans l’autre. Pour cela, la solution est d’anticiper autant que possible et de reposer au maximum sur la validité a priori. C’est-à-dire, que si on pense que l’odds ratio réel est très fort, alors on aura anticipé ce problème et on calculera le nombre de sujets nécessaires afin d’avoir suffisamment d’exposés dans les deux groupes.

Approche combinée

L’approche que je conseille combine et hiérarchise les deux méthodes : a priori et a posteriori. Idéalement, il faut anticiper au maximum la validité statistique, en se mettant de bonnes marges de sécurité. On n’est néanmoins jamais à l’abri de résultats très écartés de ce qu’on avait anticipé. La vérification a posteriori de la validité est avant tout un filet de sécurité : sa présence est importante, mais il doit servir le plus rarement possible. C’est ce que les anglo-saxons appelleraient un sanity check. En cas de chute dans la filet de sécurité, je ne donnerai pas de règle quant aux analyses à réaliser. Une ou plusieurs choses se sont complètement écartées de ce qu’on anticipait. D’un côté, toute analyse que l’on réalisera sera biaisée puisque le schéma fréquentiste global, supposant que la procédure de l’analyse est parfaitement prédéfinie, sera violé. D’un autre côté, c’est aussi une opportunité d’avancer vraiment dans la connaissance sur le sujet. Les échecs sont au moins aussi importants que les succès en recherche. Le processus de publication, obligeant tous les articles à être les plus lisses et prévisibles possibles ne facilite malheureusement pas la publication de ces résultats pourtant nécessaires à l’avancée de la recherche.

On remarquera que l’anticipation, bien que parfois délicate, facilite beaucoup l’identification des situations où les résultats sont écartés de ce qu’on avait anticipé. Cela limite le HARKing (Hypothesis After Results are Known) et permet de vraiment apprendre des résultats.

Méthode générale a posteriori

Les dysfonctionnements pouvant conduire à une invalidité statistique détectable a posteriori sont:

  1. Hypothèses erronées sur la population : prévalences, incidences, distributions ou paramètres fortement différents de ce qui avait été anticipé
  2. Mauvaise connaissance de l’estimateur, qui, même dans des conditions a priori bonnes, s’avère être biaisé

Il existe une méthode très générale de validation a posteriori qui permet de limiter beaucoup ce second risque : le bootstrap.

Par exemple, pour évaluer la validité de l’intervalle de confiance de Student, on crée des rééchantillonnages de l’échantillon de taille N, par tirage au sort de N observations avec remise parmi les N observations. On calcule l’intervalle de confiance de Student dans chacun des N rééchantillonnages puis on estime la proportion des intervalles de confiance entièrement en dessous et au-dessus de la différence de moyennes observée sur l’échantillon initial. On compare ces proportions empiriques aux proportions nominales : généralement 2,5% de risque de sous-estimation et 2,5% de risque de surestimation.

Cette méthode est aussi adaptée aux estimateurs plus complexes, par exemple, dans des situations où un score de propension ou un modèle multivarié est impliqué. Dans ces situations il est important de bootstrapper la procédure dans son ensemble : par exemple, pour un ajustement sur le score de propension, il faut recalculer le score de propension sur chaque échantillon.

Un bénéfice secondaire majeur de cette stratégie est la progressive est l’acquisition de la connaissance du comportement des estimateurs que vous acquerrez progresivement durant votre pratique statistique. Cela facilitera l’anticipation de la validité. Il est aussi possible de faire des variantes, par exemple, en changeant les tailles des rééchantillonnages tout en rééchantillonnant toujours depuis le même échantillon initial. Cela vous permettra d’identifier des « limites » de validité. Il est aussi possible de déformer les distributions.

C’est ainsi que vous verrez que les tests de normalité n’ont aucune pertinence pour la méthode de Student : (1) si les deux groupes sont de taille peu différentes, la méthode est particulièrement robuste (2) la méthode est très résistante aux distributions discrètes, platykurtiques ou même leptokurtiques tant que ces distributions restent symétriques et bornées (3) en présence d’outliers, c’est cette fréquence d’outlier qui est importante à connaître : on peut considérer que la validité devient superposable à celle d’un estimateur de différences de pourcentages d’outliers. Nécessité d’une quarantaine d’outliers s’ils sont tous situés dans le même groupe, mais d’un bien plus petit nombre si les deux groupes sont de taille proche et la fréquence est peu différente entre les deux groupes.

Il est aussi possible, pour sa culture personnelle, de comparer plusieurs estimateurs afin de mieux les connaître.

Limites de la méthode

Cette méthode d’estimation empirique de la validité de l’estimateur par rééchantillonnage souffre d’une limite principale : le rééchantillonnage est réalisé à partir de l’échantillon plutôt que de la population. Si on est largement dans les conditions de validité, l’échantillon devrait être suffisament grand pour que les différences population-échantillon n’aient pas une influence majeure sur les résultats, et confirment alors la validité. Si on est largement en dehors des conditions de validité, l’échantillon pourra différer de la population, mais il est probable que des problèmes de validité restent manifestement visibles par les rééchantillonnages, même s’il existe des exceptions dont nous discuterons. Le problème surviendra généralement lorsqu’on est proches des limites de validité. On pourra alors surestimer ou sous-estimer la validité de l’estimateur selon que l’échantillon est malchanceux ou chanceux. Si on applique le principe « better safe than sorry« , il faut donc se mettre une certaine marge de sécurité, en ayant un niveau d’exigence assez élevé sur l’inflation ou la déflation des risques de sous-estimation et surestimation par intervalles de confiance. On peut aussi artificiellement baisser la taille des rééchantillonnages pour se mettre en situation défavorable.

Nous allons maintenant décrire les situations d’incapacité totale de la méthode du bootstrap à identifier un problème majeur de validité. La première situation est celle des outliers cachés. En reprenant l’exemple de l’assurance des satellites, si la base de données ne contient aucune défaillance de lanceur, on n’observera aucun problème de validité et on conclura, à tort, que la méthode de Student est tout à fait valide. C’est pourquoi, ce ne sera que par connaissance a priori du problème qu’on saura qu’on est largement en dehors des conditions de validité. Par contre, si on observe un ou quelques outliers, on devrait prendre conscience du problème même si on aura du mal à l’évaluer précisément. En effet, considérant l’indépendance entre les observations, si on observe trois outliers, on pourra considérer que l’espérance du nombre d’outliers (paramètre lambda d’une loi de Poisson) est compris entre 0,24 et 7,2 (intervalle de confiance de Garwood, cas limite d’un Clopper-Pearson avec dénominateur infini). Dans les deux scenarii extrêmes, on reste en dehors des conditions de validité.

La seconde situation est la présence de discontinuités dans l’estimateur. Par exemple, un estimateur d’effet avec des ajustements réalisés ou pas selon leur degré de significativité (p.e. inclusion des variables pas à pas) aura des fluctuations d’échantillonnage complètement différentes sur des rééchantillonnages que si on répétait l’expérience depuis la population, en raison de la non reproductibilité des petits p. Par exemple, si la puissance statistique pour une covariable est à 30%, et que cette covariable est fortement liée à l’exposition d’intérêt, il est possible que les fluctuations soient cahotiques et les approximations normales complètement invalides. Sur l’échantillon, cette covariable a 20% de risque de faire un petit p bilatéral supérieur ou égal à 0.78 (calcul : pnorm(-(qnorm(0.975)+qnorm(0.20)+ qnorm(0.20)))*2) ce qui conduirait à une fréquence d’inclusion de 5,9% (calcul : pnorm(qnorm(0.20)+ qnorm(0.20)) + 1-pnorm(qnorm(0.20)+ qnorm(0.20)+2*qnorm(0.975))) ce qui va très fortement sous-estimer les fluctuations discrètes engendrées par l’inclusion aléatoire de cette covariable dans le modèle.

De toute façon, il est préférable d’utiliser des estimateurs ne présentant que des fluctuations continues, mais ce n’est pas toujours évitable à 100%, notamment lorsqu’on manipule des lois binomiales ou de Poisson dans des cas limites.

La réalité des choses : pressions extérieures

Au final, on peut être obligé de calculer des statistiques sur des résultats secondaires, en sachant, dès l’élaboration du protocole, qu’on fonce droit dans le mur, parce que les investigateurs ou évaluateurs du projet insistent sur la nécessité d’utiliser un critère de jugement. Par exemple, on m’a imposé de comparer un incidence ratio à la valeur 1, sachant, qu’en étant optimiste, il était tout au plus à 1,20. Le nombre d’événements attendu, les deux groupes réunis, était inférieur à 2. Il faudrait des centaines d’événements si on voulait une puissance acceptable. Ma stratégie, dans ces situations, est de planifier l’échec, en imposant un critère statistique conditionnant le calcul à une statistique sur l’échantillon. Par exemple, je rédigerai « en dessous de 20 événements au total, aucune comparaison ne sera faite. » dans un contexte ou je sais très bien qu’on n’atteindra pas les 20 événements. Si j’avais tort, cela me placerait dans la même situation qu’un échec inattendu : celle d’un résultat extraordinaire qui mérite publication.

Apprentissage

J’ai déjà mentionné qu’en évaluant a posteriori la validité des estimateurs on apprenait à les connaître. J’insiste là-dessus. De mon point de vue un bon statisticien connaît les outils qu’il utilise. La formation continue, à partir de ses propres analyses lui permet d’acquérir cette connaissance, concernant à la fois les estimateurs, mais aussi les distributions univariées et multivariées du domaine qu’il analyse. C’est ainsi qu’il pourra anticiper un maximum de choses lors de l’élaboration du protocole. Contrairement à la plupart des formations initiales et continues qui restent trop générales, cette formation sur le tas permet l’acquisition d’une quantification fine et intuitive des distributions les plus courantes ainsi que du comportements des estimateurs.

Erreur fréquente : ajustement sur baseline

Un billet pour vous présenter une erreur d’analyse tellement fréquente que c’est l’analyse correcte qui devient une exception. C’est le cas de la comparaison de scores quantitatifs (p.e. scores de qualité de vie, scores fonctionnels) après une intervention (p.e. 3 mois après) entre deux groupes de sujets susceptibles de différer sur le score avant l’intervention (baseline). Le cas le plus problématique est celui où les deux groupes ne sont pas randomisés. Cela peut correspondre à l’analyse de facteurs pronostiques sous traitement, les groupes étant définis par ces facteurs pronostiques, ou à la comparaison non randomisée de deux interventions.

Même si une partie des réflexions s’appliquent aussi bien à la recherche de facteurs pronostiques sous traitement qu’à la comparaison thérapeutique, je me focaliserai sur ce second scenario car il est difficile de définir si un pronostic (évolution favorable ou défavorable) est meilleur qu’un autre lorsque le point de départ (état initial) diffère. Pour les facteurs pronostiques, les différences à baseline ne sont pas forcément accidentelles mais sont souvent réelles, devant être prises en compte pour l’interprétation, alors que pour les comparaisons thérapeutiques, ces différences sont typiquement dues à des biais d’indication et pourraient être rectifiées par un essai clinique randomisé.

Les trois manières les plus fréquentes de comparer les deux groupes sur leur score post-intervention sont:

  1. La comparaison brute des scores post-intervention, sans ajustement (p.e. test de Student sur séries indépendantes)
  2. La comparaison des changements (aussi appelées deltas ou différences appariées) de score entre les deux groupes. Ça peut être fait par un test de Student sur séries indépendantes sur les deltas, ou par un test d’interaction temps×groupe dans un modèle linéaire à effets mixtes.
  3. Comparaison des scores post-intervention entre les deux groupes avec ajustement sur le score pré-intervention (modèle linéaire général).

La méthode N°3 est la seule a être correcte dans les études comparant deux thérapeutiques. La méthode N°1 ignore les différences pré-intervention et est donc sujette aux facteurs de confusion tels que le biais d’indication qui tendrait, par exemple, à faire donner le traitement innovant aux patients dont le pronostic est meilleur, ou pire, que les autres. La méthode N°2 donne l’illusion de corriger les différences pré-existantes au temps pré-intervention mais en réalité, sur-corrige systématiquement la différence. Le groupe avantagé à baseline se trouvera alors fortement désavantagé. Un biais d’indication positif deviendra un biais d’indication négatif.

Pour comprendre ça, il faut comprendre comment fonctionne la méthode N°3, pourquoi elle est correcte, et comment les méthodes N°1 et N°2 se comportent par rapport à cette méthode de référence.

Régression vers la moyenne

Il faut comprendre que même pour une maladie chronique « stable », deux mesures répétées ne sont jamais parfaitement corrélées, ne serait-ce qu’à cause des erreurs de mesure et fluctuations d’échantillonnage de chacune des deux mesures. Si la maladie est stable, sans la moindre intervention, alors une régression linéaire d’une mesure faite à T1 sur une mesure faite à T2 (p.e. trois mois plus tard), va trouver une pente systématiquement inférieure à 1.

Ce petit exemple de code R illustre mon propos:

vraie_valeur_T1=rnorm(100000) # valeur réelle de l'état du sujet
vraie_valeur_T2=vraie_valeur_T1

mesure_T1 = vraie_valeur_T1 + rnorm(100000)
mesure_T2 = vraie_valeur_T2 + rnorm(100000)

lm(mesure_T2 ~ mesure_T1) # pente à 0.50

Dans l’exemple ci-dessus, je distingue deux concepts:

  1. la valeur du score du sujet, qui reflète son état général de base, tel qu’on pourrait l’obtenir en faisant la moyenne d’un grand nombre de mesures étalées sur une période de temps suffisante
  2. la mesure du score du sujet, qui reflète la valeur du score à laquelle s’ajoutent la fluctuation d’échantillonnage attribuable au jour de la mesure ainsi que l’erreur de mesure. Ces sources de variation sont imprévisibles et seront donc appelées l’aléatoire.

Dans l’exemple ci-dessus, bien que les valeurs à T1 et T2 soient identiques pour chacun des sujets, la pente de la régression linéaire faisant le lien entre les deux mesures est à 0.50.

Comment interpréter cette pente de régression à 0.50 ? Cela veut dire que l’espérance de la mesure à T2, conditionnelle à une mesure à T1, est égale à la moyenne générale plus 0.50 fois la mesure à T1. C’est-à-dire, un sujet dont la mesure à T1 est à +1 écart-type de la moyenne générale, devrait se trouver, en moyenne à +0.5 écart-type de la moyenne générale à T2. Un sujet qui se trouve à -1 écart-type de la moyenne générale à T1, devrait se trouver, en moyenne à -0.5 écart-type de la moyenne générale. Toute mesure écartée de la moyenne générale à T1, devrait se trouver, en moyenne, plus proche de la moyenne générale à T1. Ce phénomène est connu sous le nom de « régression vers la moyenne ». L’explication intuitive de ce phénomène, c’est qu’une mesure à +1 écart-type à T1 est explicable, en partie par le fait que le sujet à une valeur à T1 plus grande que les autres (ici, +0.5 écart-type en moyenne), mais aussi parce qu’il a eu une source de variation aléatoire plus grande que les autres (ici +0.5 écart-type en moyenne). Comme l’aléatoire est, par principe, indépendant de tout ce qui est observable, lors de la deuxième mesure à T2, l’aléatoire sera en moyenne nul (légèrement positif ou négatif sur une observation particulière), de telle sorte qu’en moyenne, il ne restera que les +0.5 écarts-types attribuables à une valeur réellement plus élevée à T1, alors qu’en moyenne l’aléatoire sera revenu à 0.

NB : Cette idée qu’une valeur à T1 à +1 écart-type de la moyenne générale est attribuable en partie à une différence réelle de valeur, et en partie attribuable à de l’aléatoire, est vraie, en moyenne seulement. C’est-à-dire, qu’à un niveau individuel, il est possible qu’une mesure à +1 écart-type soit la somme d’une valeur réelle à -0.2 écart-type et d’un aléatoire à +1.2 écart-type. Néanmoins, en moyenne, le phénomène décrit ci-dessus s’appliquera.

L’exemple fourni est assez artificiel, avec une pente à 0.50, explicable par le fait que la variance inter-sujet de la valeur réelle est égal à la variance de l’aléatoire. Si la part d’aléatoire est plus importante, alors la pente va diminuer, conduisant à une régression vers la moyenne plus violente. Au contraire, si la part d’aléatoire est minime, on tendra vers une pente à 1.

ratio = 3 # ratio entre l'écart-type aléatoire et l'écart-type
# de la vraie valeur

vraie_valeur_T1=rnorm(1000000)
vraie_valeur_T2=vraie_valeur_T1

mesure_T1 = vraie_valeur_T1 + ratio*rnorm(1000000)
mesure_T2 = vraie_valeur_T2 + ratio*rnorm(1000000)

lm(mesure_T2 ~ mesure_T1) # pente à 0.10

Dans l’exemple ci-dessus, le rapport des écarts-types est à 3, le rapport des variances est alors à 9 et la pente est à 1/(9+1) = 0.10. Le phénomène de régression vers la moyenne est très fort, parce que l’aléatoire est très fort; ainsi, une mesure très élevée témoigne très certainement d’une grosse part d’aléatoire plutôt que d’une valeur réellement très forte.

Régression vers la moyenne sous traitement

Comme décrit ci-dessus, lorsqu’il y a des mesures répétées en état stable, la deuxième mesure tendra, en moyenne, à se rapprocher de la moyenne générale par rapport à la première. Lorsque les patients bénéficient d’un traitement, alors il y aura une tendance générale à l’amélioration pour deux phénomènes :

  1. Les effets contextuels. Notamment la sélection des sujets au moment où ils vont le moins bien, et donc, dont l’aléatoire est défavorable à l’inclusion. Cela peut aussi correspondre à une maladie qui tend à la guérison spontanée.
  2. L’effet réel du traitement

Cette tendance générale s’exprimera comme la différence de moyennes entre T1 et T2. Le phénomène de régression vers la moyenne persistera après soustraction de cette tendance générale aux mesures de T2. Par ailleurs, de nouveaux phénomènes apparaîtront pour les scores de symptômes. Par exemple, si l’évolution est favorable pour la majorité des sujets, il est probable que la variance des scores soit plus faible après traitement qu’avant. En effet, pour un score symptomatique, d’autant plus élevé que les symptômes sont sévères et nombreux, les sujets auront majoritairement des scores proches de zéro après traitement. Dans ce cas, la pente de régression expliquant la mesure à T2 par la mesure à T1 sera encore plus proche de zéro. Plus rarement, lorsque l’effet du traitement est très variable, avec des sujets répondeurs et des non répondeurs, il est possible que la variance augmente, mais encore une fois, la pente de régression devrait rester inférieure à 1. D’un point de vue théorique, cette pente pourrait dépasser 1, dans des conditions qui ne sont jamais rencontrées dans la pratique : faibles erreurs de mesures et très faibles fluctuations d’échantillonnage et réponse au traitement extrêmement corrélée à la valeur à baseline, avec une variance à T2 largement supérieure à la variance à T1.

Ajustement sur la mesure pré-intervention

Comment comparer deux traitements sur une mesure post-intervention (T2) lorsque les mesures pré-intervention (T1) diffèrent ?

L’idée, c’est de rectifier les différences à T1 telles qu’elles se répercutent à T2. Cela se fait par une estimation empirique de la pente de la régression de la mesure à T2 par la mesure à T1. Cela se fait très bien dans un modèle linéaire général expliquant la mesure à T2 par l’effet groupe (tt innovant vs tt de référence) et par la mesure à T1.

Ci-dessous un exemple illustratif, dans lequel il existerait un biais d’indication majeur où les patients ayant le plus de symptômes (score élevé) bénéficieraient du traitement innovant alors que les patients ayant le moins de symptômes auraient le traitement de référence:

set.seed(2021) # graine pseudo-aléatoire pour avoir des résultats reproductibles

effet_tt = -1 # même effet traitement dans les deux groupes
N1 = 5000000 # taille du groupe 1
N2 = 5000000 # taille du groupe 2
N = N1+N2
ratio = 3 # ratio entre l'écart-type aléatoire et l'écart-type de la vraie valeur

vraie_valeur_T1=rnorm(N)
mesure_T1 = vraie_valeur_T1 + ratio*rnorm(N) # on ajoute l'erreur de mesure/fluctuation aléatoire
# le traitement innovant est donné aux sujets avec plus de s
ymptômes
tt_innovant = (mesure_T1 > 1.00) 

# peu importe le groupe, l'effet du traitement est le même dans les deux groupes
# car le traitement innovant n'est pas meilleur que le traitement de référence
vraie_valeur_T2=vraie_valeur_T1 + effet_tt

mesure_T2 = vraie_valeur_T2 + ratio*rnorm(N) # on ajoute l'erreur de mesure/fluctuation aléatoire

mean(mesure_T1[tt_innovant]) # mesure à T1 en moyenne à 3.19 dans le groupe tt innovant
mean(mesure_T1[!tt_innovant]) # mesure à T1 en moyenne à -1.92 dans le groupe tt de référence
mean(mesure_T2[tt_innovant]) # mesure à T1 en moyenne à -0.68 dans le groupe tt innovant
mean(mesure_T2[!tt_innovant]) # mesure à T1 en moyenne à -1.19 dans le groupe tt de référence

# le modèle correct
# trouve une différence nulle entre traitement innovant et traitement de référence
# et trouve un effet (pente) à 0.10 de la mesure à T1 sur la mesure à T2
lm(mesure_T2 ~ mesure_T1+tt_innovant)

Dans cet exemple, il existe un fort effet de régression vers la moyenne et quand bien même les deux traitements sont équivalents, la différence entre les groupes de traitement que l’on trouve à T1 (moyenne à 3.19 vs -1.92, différence à 5.1) se réduit par un facteur dix à T2 (-0.68 vs -1.19, différence à 0.51). Ce facteur dix est explicable par une pente de régression de la mesure à T1 sur la mesure à T2 qui est égale à un dixième (0.10), révélant une forte régression vers la moyenne. Le modèle linéaire ajusté sur la mesure à T1 est tout à fait capable de comparer les valeurs à T2 (-0.68 vs -1.19) en rectifiant la différence attribuable à la différence de moyennes à T1 (3.19 vs -1.92). Cela permet de conclure à l’équivalence des traitements (ici, l’échantillon étant immense, on peut considérer que tout est exact) car la différence à T2 n’est pas supérieure à ce qui est attendu (0.51 unité), compte tenu de la différence à baseline (5.1 unités) et de la pente de régression de T1 sur T2 (pente à 0.10).

Autres ajustements corrects ou inappropriés

Cependant, si on comparait les moyennes des mesures à T2 avec un test de Student sur séries indépendantes, en ignorant ainsi, les différences à T1, on conclurait que les sujets avec traitement innovant ont un plus grand nombre de symptômes (+0.51) que les sujets avec le traitement de référence. Ainsi, on conclurait, à tort, que le traitement innovant est moins bon que le traitement de référence alors qu’ils sont en réalité équivalents. Cet effet à +0.51 est due au biais d’indication.

Enfin, si on comparait les moyennes des changements (différences appariées T2 moins T1), on conclurait à tort que le traitement innovant est très largement supérieur au traitement de référence. En effet, cette différence de différences est égale à 0.51-5.1 = -4.6, témoignant d’une réduction des symptômes beaucoup plus marquée dans le groupe de traitement innovant bien qu’en réalité les deux traitements soient équivalents.

Comme dit en introduction, la comparaison des moyennes de changements conduit à une sur-correction des différences à baseline, transformant le désavantage du groupe traité par le traitement innovant, en un avantage.

Un autre point de vue, sur ces deux analyses biaisées, c’est de les interpréter comme une analyse ajustée dans lesquelles la pente de régression de la mesure T2 sur la mesure T1 est artificiellement imposée à zéro (pour le Student à T2 sans ajustement sur la baseline) ou artificiellement imposée à 1 (pour la comparaison de moyennes des changements). Cela peut se décrire par le code R suivant, complétant le précédent:

lm(mesure_T2 ~ offset(mesure_T1*0)+tt_innovant) # modèle ignorant les différences à baseline (Student, sans ajustemnet)
lm(mesure_T2 ~ offset(mesure_T1*1)+tt_innovant) # pente imposée à 1, comparant les moyennes des changements

En effet, les deux équations suivantes, de modèles linéaires, sont équivalentes:

T2 = beta0 + beta1×tt + 1×T1 + epsilon

(T2-T1) = beta0 + beta1×tt + epsilon

Où epsilon suit une loi normale de moyenne zéro et écart-type sigma.

On remarquera aussi que l’analyse des changements, ajustée sur les mesures à T1, est équivalente à l’analyse des mesures à T2, ajustée sur les mesures à T1. C’est-à-dire qu’on obtient exactement le même effet traitement avec les deux modèles suivants:

lm(mesure_T2 ~ mesure_T1+tt_innovant)
lm(mesure_T2-mesure_T1 ~ mesure_T1+tt_innovant)

Encore une fois, ça découle assez directement des équations des modèles linéaires:

T2 = beta0 + beta1×tt + beta2×T1 + epsilon

(T2-T1) = beta0′ + beta1’×tt + beta2’×T1 + epsilon

On constate que les deux équations sont équivalentes pour beta0=beta0′, beta1=beta1′ et beta2’=beta2-1. En bref, les coefficients de l’ordonnée à l’origine (intercept) et de l’effet traitement sont inchangés; seul change le coefficient de la pente, qui est réduit d’une unité pour la seconde analyse par rapport à la première. Autrement, les modèles sont équivalents, au sens où ils prédisent exactement les mêmes valeurs pour les mêmes sujets. Les résidus des modèles sont identiques. Les coefficients beta0 et beta1 ont exactement les mêmes estimations ponctuelles sur n’importe quel échantillon.

C’est très intéressant parce que ça permet de présenter exactement les mêmes résultats sous un angle plus plaisant pour les investigateurs qui ne comprennent rien aux phénomènes de régression vers la moyenne et qui veulent absolument comparer les moyennes des changements : on peut leur dire qu’on a comparé les moyennes des changements, comme ils le souhaitaient, mais qu’on a ajusté sur la mesure à T1.

Modèle à effets mixtes

Comme présenté en introduction, le modèle à effets mixtes expliquant à la fois la mesure à T1 et à T2, avec un intercept aléatoire sujet, un effet temps, un effet groupe et une interaction temps×groupe, est équivalent à une comparaison de moyennes des changements (différences appariées) et conduit donc a une sur-correction de la différence à baseline. En l’absence de données manquantes, les estimations ponctuelles sont identiques. Néanmoins, ces modèles à effets mixtes rajoutent un grand nombre d’hypothèses dont on doit réfléchir aux conséquences. D’abord, il y a une hypothèse d’homoscédasticité (égalité des variances résiduelles) à T1 et à T2 alors que le simple Student fait seulement l’hypothèse d’homoscédasticité entre les groupes sur les différences appariées. Heureusement, comme le nombre de mesures à T1 et à T2 sont égales, en l’absence de données manquantes, les écarts à cette hypothèse de validité semblent n’engendrer aucun biais sur les intervalles de confiance. Il y a aussi une hypothèse de normalité des résidus qui semble, comme le Student, être relâchée lorsque l’échantillon est suffisamment grand. Le conditionnement à l’exposition, notamment à l’effet sujet, ne pose pas non plus de problème lorsque le modèle est linéaire avec fonction de lien identité. On peut aussi redouter que la variabilité de l’effet du traitement d’un sujet à l’autre pose problème, mais en réalité, cela est intégré dans le résidu, et se contente d’augmenter la variance résiduelle à T2, avec des problèmes d’hétéroscédasticité qui s’avèrent n’avoir pas de conséquence. En l’absence de données manquantes, c’est donc, équivalent, en estimation ponctuelle, et en incertitude, à un Student sur les différences appariées. En pratique, il y a quelques minimes biais supplémentaires, tels que le décompte des degrés de libertés qui est souvent incorrect dans les modèles à effets mixtes (dépend du logiciel et du package utilisé). En présence de données manquantes, le modèle à effets mixtes donne l’illusion de faire participer les perdus de vue, alors qu’en réalité, ils ne participent quasiment pas à la statistique car ils ne fournissent qu’une information infime sur la différence appariée ; et cette information repose énormément sur l’hypothèse de normalité. Au mieux, ils servent un tout petit peu à estimer la variance, en risquant de la biaiser un tout petit peu ; mais même là, leur variance étant la somme de la variance inter-sujet et intra-sujet, ils ne fournissent quasiment aucune information sur l’un ou l’autre des deux ; et cette information repose aussi massivement sur l’hypothèse de normalité. En bref, c’est une manière d’exclure les données manquantes sans le dire.

Essais cliniques randomisés

Pour les essais cliniques randomisés bien menés, la différence à baseline a une espérance nulle et les trois analyses ont la même espérance : comparaison sans ajustement à T2, comparaison ajustée sur T1 et comparaison des différences T2-T1. Néanmoins, la variance résiduelle sera minimale avec la comparaison ajustée sur T1, la précision statistique sera maximale et la puissance sera maximale, sauf si la pente de régression de T1 sur T2 est extrêmement proche de 0 ou de 1. J’affectionne particulièrement la comparaison à T2 sans ajustement, car cela est l’analyse la plus simple possible. Cela évite d’être accusé de p-Hacking.

D’ailleurs, c’est un des problèmes de ces multiples possibilités d’analyse : il est très facile de p-hacker. Même si les analyses sont équivalentes en espérance, sur un échantillon donné, elles peuvent différer de manière substantielle, fournissant alors une grande liberté de choix du petit p.

Erreur fréquente : facteurs de réponse au traitement

Un petit mot pour dire qu’il ne faut pas confondre pronostic sous traitement et réponse au traitement. La réponse à un traitement, c’est la différence d’évolution entre un sujet qui aurait le traitement et le même sujet qui ne l’aurait pas. Généralement, on peut évaluer la réponse moyenne à un traitement par un essai clinique randomisé en groupes parallèles. La réponse individuelle est bien plus difficile à évaluer, voire impossible, car il n’est pas forcément possible de savoir qu’elle aurait été l’évolution si le patient n’avait pas reçu la prise en charge qu’il eût.

Ce qu’il faut faire

L’identification d’un facteur de réponse au traitement nécessite schématiquement quatre groupes:

  1. un groupe exposé au facteur et prenant le traitement
  2. un groupe exposé au facteur et ne prenant pas le traitement
  3. un groupe non exposé au facteur et prenant le traitement
  4. un groupe non exposé au facteur et ne prenant pas le traitement

(pour un facteur quantitatif, on peut distinguer différents niveaux d’exposition, mais l’idée reste la même, il doit y avoir une variance de ce facteur d’exposition dans un groupe traité comme dans un groupe non traité)

La différence moyenne entre (1) et (2) représente la réponse moyenne des exposés alors que la différence moyenne entre (3) et (4) représente la réponse moyenne des non exposés. La différence des différences (interaction) représente la différence de réponses et si elle est statistiquement significativement différente de zéro (ou d’un seuil de significativité clinique) conduit à la conclusion que l’exposition considérée est un facteur de réponse au traitement.

Ces analyses de facteur de réponse au traitement correspondent aux fameuses analyses en sous-groupes que l’on retrouve dans bon nombre d’essais cliniques randomisés.

Ce qu’il ne faut pas faire

Malheureusement la majorité des soi-disant recherche de facteurs de réponse à un traitement qu’il m’ait été donné de voir ne contenaient que deux groupes: les patients traités et exposés au facteur et les patients traités et non exposés au facteur. Cela permet seulement de trouver des facteurs pronostics sous traitement, c’est-à-dire des facteurs prédictifs d’une évolution favorable ou défavorable, et ce, sur une population de patients tous traités.

Comme exemple frappant, je suggère de considérer l’analyse de soi-disant facteurs de réponse à l’homéopathie sur le cancer de la prostate. Le critère de jugement principal serait la survie globale de patients. On constaterait alors que les patients peu symptomatiques, dont le stade tumoral est précoce (notamment sans métastases) et de bas grade Gleason sont les meilleurs répondeurs à l’homéopathie.

Régression vers la moyenne

Pourquoi confond-on la réponse au traitement et le pronostic sous traitement ? Pourquoi toute évolution clinique n’est-elle pas totalement attribuable au traitement ?

Certaines maladies ont une évolution spontanée tendant vers une évolution inexorablement défavorable; cela comprend notamment la plupart des maladies neurodégénératives telles que les maladies démentielles, la maladie de Parkinson, la sclérose latérale amyotrophique. Même si cette évolution défavorable ne concerne pas forcément 100% des patients, on peut dire que l’état clinique moyen se dégrade.

D’autres maladies ont une évolution moyenne plus stable, comme l’asthme chronique. Néanmoins, il est très rare que la stabilité individuelle soit parfaite. Il existera presque toujours des fluctuations intra-sujet, c’est-à-dire, des jours, semaines ou mois plus symptomatiques que d’autres. Le phénomène de régression vers la moyenne s’applique alors à ces sujets. En moyenne un sujet sélectionné sur un état clinique moins bon que la moyenne, tendra à en voir une amélioration dans les jours, semaines ou mois avenir. Or, l’initiation des traitements est rarement faite au moment où le sujet va le mieux; au contraire, ce sera lors d’un état clinique médiocre ou d’une dégradation de l’état habituel. Parfois cela est même formalisé explicitement dans les critères d’inclusions, sous forme d’un seuil sur un score ! Par simple régression vers la moyenne, l’état clinique du sujet tendra à s’améliorer quel que soit le traitement donné; c’est une des principales raisons justifiant la perception d’efficacité de l’homéopathie. Il s’agit alors de l’évolution spontanée de la maladie.

Si on analyse la différence entre deux mesures répétées d’un score des symptômes, à baseline et un certain temps après (p.e. 3 mois), alors il s’agit de la somme de deux composantes:

  1. L’évolution spontanée
  2. La réponse au traitement

Malheureusement, c’est souvent interprété comme la réponse au traitement.

Il apparaît alors que les patients dont l’état est le pire vont le mieux s’améliorer, ce qui pourra être, à tort, considéré comme une meilleure réponse au traitement.

Si on a en conscience cela, il faut savoir qu’on peut avoir une réponse au traitement quand bien même l’état clinique du patient est stable, comme pour la sclérose latérale amyotrophique pour laquelle l’évolution naturelle est très défavorable et une stabilisation prolongée déjà acceptable pour un traitement. Et bien sûr, il peut y avoir une évolution spontanée favorable sans la moindre réponse au traitement. Un traitement homéopathique d’une angine tendra à une résolution totale des symptômes en généralement quelques jours, ainsi que l’absence totale de traitement. La réponse doit être évaluée par la différence entre les des deux effets.

C’est peut-être pour cette confusion entre réponse au traitement et pronostic sous traitement que l’homéopathie est aussi populaire en France.

Pour aller plus loin

Comme mentionné au début, un facteur de réponse au traitement est un facteur ayant une interaction statistique avec l’effet du traitement, c’est-à-dire, la différence entre un groupe traité et un groupe non traité. La notion d’interaction est assez simple lorsque les états cliniques sans traitement sont identiques chez les exposés et les non exposés au facteur. Par contre, si ces états diffèrent, l’interaction statistique dépend alors du modèle employé, pouvant conduire à des conclusions opposées.

Considérons, par exemple, que l’on s’intéresse à un critère de jugement binaire, tel que « mortalité à 1 an ».

Non traitésTraités
Non exposés30%10%
Exposés60%30%
Taux de mortalité à 1 an selon la présence d’exposition et l’usage d’un traitement spécifique

Pour une différence absolue de mortalité, on peut dire que les non exposés ont une meilleure réponse au traitement (-30% de décès) que les non exposés (-20%). Pour un risque relatif de décès, au contraire, les non exposés ont une réponse moindre (RR = 0.50) aux exposés (RR = 0.33). Pour un odds ratio de risque de décès, c’est à nouveau les non exposés qui ont une meilleure réponse (0.26 vs 0.29). La réponse pourrait encore différer si l’on s’intéressait aux rapport des chances de survie plutôt que des risques de décès. Bien sûr, les hazard ratio conduisent encore à des conclusions différentes.

Est-ce vraiment si important que ça ? De toute façon, les deux groupes, exposés et non exposés, bénéficient chacun du traitement et ces deux groupes ne sont pas comparables de toute façon. Au final, à quoi sert-il de comparer leur réponse ? Il pourrait finalement être plus sage de juste comparer leur pronostic sous traitement ou de ne rien comparer du tout.

Là où il est crucial d’identifier des facteurs de réponse au traitement, c’est lorsque certains sous-groupes ne répondent pas du tout, voire pire, ont une réponse négative au traitement (interaction qualitative). Dans ces situations, la statistique d’analyse de l’effet du traitement n’importe pas, puisque toutes tendront à la même conclusion.

On peut aussi s’interroger sur la pratique consistant à tenter de prouver l’homogénéité de l’effet d’un traitement dans divers sous-groupes en formulant une hypothèse nulle d’absence d’interaction, puis en acceptant cette hypothèse après réalisation d’un test statistique sous-puissant. Comme précédemment décrit, les interactions quantitatives, consistant à des effets plus ou moins forts du traitement selon les sous-groupes, sont bien moins préoccupantes que les interactions qualitatives; mieux vaudrait juste prouver l’existence d’un effet positif du traitement dans chacun des sous-groupes plutôt que de comparer les effets entre eux ; mais cette méthode mettrait en évidence l’incapacité à conclure sur certains sous-groupes en situation de sous-puissance statistique plutôt que de fournir une confortable acceptation de l’hypothèse nulle en situtation de sous-puissance majorée.

Sur ce, je vous laisse méditer…

Choix du critère de jugement

Incorporation de l’exposition dans le critère de jugement

À la troisième observation d’une même erreur méthodologique, dans mon expérience personnelle, je considère qu’un problème existe. Cela ne permet pas d’en évaluer la fréquence, mais cela veut dire qu’il vaut peut-être la peine de la mentionner.

Un critère de jugement principal ne doit pas être confondu avec l’exposition. Considérons une maladie évoluant sous forme de poussées plus ou moins récurrentes ou sous forme permanente pour laquelle on obtient assez souvent une réponse complète à moyen terme, avec disparition complète des symptômes, avec le traitement de référence. On compare dans un essai clinique randomisé en ouvert le traitement de référence A au traitement innovant B plus ou moins combiné au traitement A, selon le souhait du médecin. Considérons pour critère de jugement principal la « réponse complète off therapy » définie par l’absence de symptômes combinée à l’interruption totale du traitement de référence A. Le critère de jugement principal comprend, dans sa définition, l’exposition au traitement A ! On peut alors conclure que les sujets randomisés dans le groupe A vont prendre généralement le traitement A, alors que les sujets randomisés dans le groupe B±A vont le prendre moins souvent.

Afin de ne pas faire cette erreur, je vous conseille deux stratégies.

  1. D’abord, systématiquement vérifier que votre critère de jugement n’incorpore pas l’exposition
  2. Ensuite, considérer le scenario virtuel où on comparerait le traitement à lui même mais en version de couleur différente (comprimé rouge versus bleu, identique en principe actif), afin de vérifier que le traitement ne montrerait pas sa supériorité à lui-même

Critère de jugement différent selon le groupe

J’ai observé ce problème à la fois dans une étude non publiée et dans une étude publiée.

Étude publiée dans le JAMA

Pour commencer, Salminen et al, 2015 (https://dx.doi.org/10.1001/jama.2015.6154) présente un essai clinique randomisé comparant la chirurgie à l’antibiothérapie dans le traitement de l’appendicite aiguë non compliquée. Dans le groupe chirurgie, le succès (critère de jugement principal) est défini par le fait que le patient a été opéré (taux attendu de succès ~= 100%) alors que dans le groupe antibiotiques il est défini par le fait que le patient n’a pas été opéré ! La non-infériorité devait être démontrée avec une marge de -24% pour les antibiotiques vs chirurgie.

Le succès du groupe chirurgie est, par définition, un échec du groupe antibiotiques et vice versa ! Si on comparait la chirurgie à elle même, on arriverait aisément à la conclusion que la chirurgie est très largement inférieure à elle même (0% vs 100% de succès) tout en étant très largement supérieure à elle même (100% vs 0% de succès).

On peut d’ailleurs se demander comment on peut ne pas avoir 100% de succès de la chirurgie ! Outre les décès per-opératoires (très rares pour une appendicite aiguë non compliquée) il y a les annulations opératoires pour l’analyse en intention de traiter. Cela conduit à la mention « The patient randomized to appendectomy who did not have an operation had resolution of symptoms before the operation could be performed ». Peut-être était-ce juste un patient bien constipé pour lequel on a fait une erreur diagnostique ? Je ne suis pas sûr que le patient considère lui-même que c’est une évolution défavorable de sa maladie, par contre, c’est peut être une évolution défavorable des finances du chirurgien s’il travaille dans le secteur privé.

Étude non publiée

Plusieurs stratégies thérapeutiques étaient analysées, dans une optique de désescalade thérapeutique de la chimiothérapie anti-cancéreuse. Dans un sous-groupe bien spécifique les patients étaient randomisés en un groupe d’abstention thérapeutique alors que l’autre groupe bénéficiait d’une chimiothérapie. Le critère de jugement principal était la survie sans rechute dans le groupe avec chimiothérapie et de la survie sans re-progression après rechute dans le groupe d’abstention thérapeutique. C’est-à-dire que dans le groupe chimiothérapie, on mesure le délai avant une première rechute, alors que dans le groupe d’abstention thérapeutique on attend la première rechute avant de mesurer le délai entre la première et la seconde rechute.

Encore une fois, si on comparait la chimiothérapie à elle même, on trouverait une différence, parce qu’il n’y a pas de raison que le délai avant première rechute soit identique au délai entre première et seconde rechute.

Pondération d’un modèle à effets mixtes

Les modèles à effets mixtes sont très largement utilisés, mais la statistique qu’ils estiment n’est pas forcément correctement interprétée. Alors que tous les statisticiens connaissent l’interprétation de la médiane ou la moyenne, on parlera plutôt d’effet dans un modèle à effets mixtes sans forcément savoir ce que ça représente.

Pour aider à l’interprétation et au choix des statistiques, je vous propose de partir de deux exemples de mesures répétées très simples.

1er exemple : paires de jumeaux

Supposons qu’on souhaite évaluer les conséquences de la grande prématurité (< 32 semaines d’aménorrhée) sur le développement psychomoteur de l’enfant. On comparerait alors, certains outcomes de développement psychomoteur, entre enfants prématurés et non prématurés.

Plusieurs problèmes statistiques apparaissent:

  1. L’exposition (prématurité ou non) des jumeaux est parfaitement corrélée (R=1) puisque le terme de naissance est la même pour les deux jumeaux d’une paire (sauf rare exception)
  2. Le devenir des jumeaux est aussi corrélé, car ils partagent très fortement leur environnement pré-natal et post-natal et partagent plus ou moins fortement leur génôme
  3. Il existe une corrélation forte entre la gémellarité et la prématurité; c’est-à-dire, qu’un plus grand nombre de jumeaux sera retrouvé dans le groupe prématuré
  4. À prématurité égale, le devenir psychomoteur moyen d’un enfant pourrait différer selon que la grossesse soit gémellaire ou pas. Une prématurité à 31 SA chez des jumeaux peut être principalement attribuable à la gémellarité alors qu’elle sera plus souvent due à d’autres comorbidités (p.e. macrosomie, RCIU, anomalies congénitales) dans des grossesses simples. Ainsi, il existe potentiellement une corrélation entre la taille du cluster et son devenir.

Néanmoins, avant de nous concentrer sur tous ces problèmes, considérons que nous souhaitons répondre à la question : quelle est la valeur moyenne de l’échelle de développement psychomoteur à l’âge corrigé de deux ans chez les prématurés ?

Pour répondre à cette question, on doit se demander : est-ce que la moyenne porte sur les enfants, avec une paire de jumeaux qui compte pour deux enfants, ou est-ce que la moyenne porte sur les grossesses, avec une paire de jumeaux qui ne compte que pour un ? Dans le premier cas, on ferait la moyenne brute des résultats de tous les enfants, sans prendre en compte le fait que deux enfants peuvent appartenir à une paire de jumeaux. Dans le second cas, on commencerait par faire la moyenne dans chacune des paires de jumeaux, afin de n’obtenir qu’une seule valeur par paire, avant de faire la moyenne de tous les clusters. On peut aussi considérer que, dans le premier cas, on fait la moyenne non pondérée de tous les résultats de tous les enfants, alors que dans le second cas, on fait la moyenne pondérée par l’inverse du nombre d’enfants dans la grossesse.

On peut donc résumer la question à : est-ce que les deux jumeaux d’une même grossesse comptent pour deux ou ne comptent que pour un ?

Même si la réponse à cette question peut sembler difficile au premier abord, elle ne fait pas de doute pour moi. Que des séquelles graves soient retrouvées chez deux singletons issus de deux grossesses différentes ou alors qu’elles soient retrouvées chez deux jumeaux issus d’une même grossesse, le handicap populationnel conséquent est le même. Les deux jumeaux sont deux individus distincts dont la valeur de la vie compte tout autant que celle de deux individus différents. La pondération par un-demi des jumeaux supposerait que la valeur de la vie des jumeaux est deux fois moins importante que celle des singletons. La réponse, pour moi est claire: deux jumeaux comptent pour deux.

Cette distinction entre les deux moyennes revêt une grande importance à cause des problèmes N°3 et N°4 évoqués ci-dessus. En effet, à cause de ces problèmes, l’espérance d’une moyenne ou d’une diffférence de moyennes diffère selon le choix de la pondération.

2ème exemple : calendrier de symptômes

Je pars d’un exemple réel. Un essai clinique randomisé sur un traitement de l’incontinence fécale. L’objectif étant d’obtenir une diminution de la fréquence des épisodes d’incontinence de selles et d’impériosités, le critère de jugement principal est basé sur le remplissage d’un calendrier, rempli quotidiennement sur 21 jours consécutifs après trois mois de traitement. Avec 21 mesures par sujet, on atténue la variance intra-sujet. On peut raisonnablement supposer que sur la période de mesure, l’état des patients sera en moyenne stable. Néanmoins, on peut craindre que certains patients ne remplissent qu’à moitié le calendrier (pe. les 10 premiers jours de la période d’évaluation), parce que le remplissage est répétitif et ennuyeux. On peut aussi craindre que ce remplissage partiel soit corrélé à la fréquence et la sévérité des symptômes.

On retrouve donc jusqu’à 21 données par patient, corrélées les unes avec les autres. La même question que pour le premier exemple se pose : un patient qui a rempli 7 jours de calendrier doit-il compter trois fois moins qu’un patient qui a rempli correctement les 21 jours ? De mon point de vue, la réponse, cette fois-ci est inversée par rapport au premier cas. Il n’y a pas de raison de donner un plus grand poids aux patients ayant entièrement complété le questionnaire qu’aux autres. On peut craindre que leur donner un plus grand poids biaise les résultats si le taux de remplissage est corrélé au contenu du questionnaire. À l’opposé, on peut espérer qu’un remplissage de bonne qualité sur les 7 premiers jours soit déjà représentatif de la période entière, et que finalement, la corrélation entre la durée de remplissage et la valeur moyenne n’ait pas d’impact sur les résultats si on commence par calculer une unique valeur moyenne par sujet avant d’en faire la moyenne sur l’ensemble des sujets.

Synthèse des deux exemples

Nous avons vu que selon la situation, le poids donné aux mesures répétées ne devrait pas être le même. Dans le premier exemple, les clusters de deux jumeaux doivent compter double par rapport aux clusters d’un singleton. Dans le second exemple, les clusters de 21 mesures (1 seul patient) devraient compter le même poids que les clusters de 7 mesures (1 seul patient) plutôt que de compter triple. Dans le second cas, le choix est absolument critique car il est susceptible d’influencer le signe de la différence entre les deux groupes dans le cadre d’un essai clinique randomisé.

Et le modèle linéaire à effets mixtes ?

Comment le modèle linéaire à effets mixtes, avec un intercept cluster aléatoire, se comporte-t-il dans les deux exemples ci-dessus ? Pour simplifier l’exemple, considérons même que l’on ne s’intéresse qu’à calculer la moyenne d’un seul groupe avec un modèle à effets mixtes sans covariable (intercept seul).

Est-ce qu’un cluster deux fois plus grand compte double ? Ou alors, chaque cluster aurait le même poids ?

La réponse est entre les deux. Ce modèle va être interprétable comme une moyenne pondérée. Un cluster de taille deux comptera plus fortement qu’un cluter de taille un, mais ce ne sera pas le double. Ce sera une valeur intermédiaire entre 1 et 2. Cette valeur intermédiaire dépendra de la force de la corrélation intra-cluster. Si les observations d’un même cluster sont très fortement corrélées, alors le poids sera proche de 1, c’est-à-dire que le poids total d’un cluster sera presque indépendant de sa taille. À l’opposé, si les observations d’un même cluster sont très faiblement corrélées, alors le poids total d’un cluster sera proche du nombre d’observations du cluster.

Reprenons l’exemple des jumeaux (exemple 1). Le modèle à effets mixtes considèrera que la valeur de la vie de deux jumeaux est d’autant plus grande que leur devenir est divergent, faiblement corrélé. À l’opposé, si les deux jumeaux partagent leur évolution, alors leur vie ne compte que pour un. Comme si la valeur d’une vie était proportionnelle à son imprévisibilité. C’est l’aléatoire d’une vie qui en ferait la valeur ! C’est ce qu’elle a d’unique. Je vous avoue que je ne suis pas branché par ce concept de snowflake, mais c’est ce que le modèle à effets mixtes fait dans votre dos quand vous lui faites confiance.

Pour aller plus loin : quels poids donne réellement le modèle à effets mixtes

Pour pousser les choses plus loin encore, dans le premier exemple des singletons/jumeaux (clusters de taille 1 ou 2), l’estimation de moyenne du devenir des enfants prématurés par le modèle à effets mixtes est extrêmement proche de l’estimation que l’on obtient par la procédure suivante:

  1. Calculer la moyenne M1 et son erreur type dans le groupe des jumeaux prématurés (pour que l’erreur type ne soit pas biaisée, on commence par moyenner les deux jumeaux de chaque paire)
  2. Calculer la moyenne M2 et son erreur type dans le groupe des singletons prématurés
  3. Calculer la moyenne de M1 et M2 pondérée par l’inverse de la variance de M1 et de M2 (inverse du carré des erreurs types)

Sur de grands échantillons, la procédure du modèle à effets mixtes converge vers celle de la moyenne pondérée par l’inverse de la variance. J’ai fait quelques simulations pour vérifier ma théorie, et on peut, en première intention, considérer les procédures comme équivalentes tant elles convergent vite.

On comprend alors l’intérêt du modèle à effets mixtes : cette pondération par l’inverse de la variance est la procédure la plus efficace statistiquement (faible variance de l’estimateur) pour faire la moyenne de M1 et M2 sous l’hypothèse que les deux moyennes de la population µ1 et µ2 sont identiques. C’est-à-dire que si la moyenne de la mesure d’un cluster est indépendante de la taille du cluster, la procédure du modèle à effets mixtes est équivalente aux autres procédures en moyenne, tout en étant plus précise statistiquement. Mais dès qu’on s’écarte de cette hypothèse, ça perd toute sa pertinence.

Pour aller plus loin : calculs d’incertitude

Jusqu’à maintenant, je n’ai pas parlé du problème de calcul des intervalles de confiance et petits p. Je me suis concentré sur l’espérance de l’estimateur ponctuel, c’est-à-dire, la statistique qui est vraiment estimée par le modèle. Je me suis concentré là-dessus, parce que je considère que c’est la clé du choix statistique. Malheureusement, je crains que le modèle à effets mixtes soit souvent utilisé en faisant le raisonnement foireux ci-dessous:

  1. La méthode de Student repose sur l’indépendance entre les observations
  2. À cause de la corrélation des mesures, Student sous-estime la variance et fournit une inférence biaisée
  3. Les modèles à effets mixtes prennent en compte cette corrélation et donc, fournissent une inférence non biaisée
  4. Donc je vais utiliser ces modèles magiques

La faille du raisonnement, c’est que la statistique du modèle à effets mixtes diffère de la statistique de moyenne simple qu’on voulait estimer. On infère sans biais, mais sur la mauvaise cible.

Les calculs d’incertitudes, sont un détail secondaire qu’on résout généralement assez simplement. D’abord, on peut faire du boostrap sur les clusters eux-mêmes. C’est une procédure extrêmement solide, adéquate notamment lorsque les corrélations intra-clusters sont extrêmement complexes. Par exemple, le modèle à effets mixtes à intercept seul fournira des résultats biaisés dans le second exemple (calendrier des symptômes) car il ne prendra pas en compte l’auto-corrélation intra-sujet, c’est-à-dire, le fait que deux jours successifs se ressemblent plus que deux jours distants. On peut aussi utiliser un estimateur sandwich sur un modèle linéaire général pondéré. Enfin, on peut estimer les variances dans les sous-groupes (pe. clusters de taille 1 et 2) puis utiliser les formules VAR(cX)=c²VAR(X) et VAR(X+Y) = VAR(X) + VAR(Y) + 2×COV(X,Y) pour calculer la variance de la moyenne des deux sous-groupes, correctement pondérée.

Conclusion

Je pense qu’au lieu de se concentrer sur des modèles, c’est-à-dire, un ensemble d’hyptohèses sur les processus de génération des données et la forme des relations, on devrait raisonner en termes de statistique : quelle valeur synthétique représente le mieux ce qui m’intéresse. Ce n’est qu’après avoir défini cette statistique que les problématiques d’échantillon fini doivent être considérés : rectification des biais d’estimateur ponctuel et calcul des incertitudes.

Il est parfois nécessaire de reposer sur des modèles, mais les conséquences de la violation des hypothèses sous-jacentes, qui est systématique, doit être connue afin de pouvoir interpréter correctement les résultats.

Calulatrices graphiques

Un petit billet sur le marché étonnant des calculatrices graphiques en 2021. Le marché orbite autour des épreuves et concours, tels que le baccalauréat en France. Les constructeurs adaptent leurs modèles au programme et aux réglementations, tels que l’obligation du mode examen (https://calculatrice-scientifique.eu/mode-examen-concours/) ou l’ajout du langage Python dans les derniers modèles. Le marché est très large, puisque rien qu’en France, environ 750 000 candidats passent le baccalauréat chaque année. Même s’il existe des spécificités nationales, les constructeurs utilisent les mêmes modèles partout dans le monde. Casio, Texas Instruments et Helwet Packard sont les principaux constructeurs.

Le prix de ce matériel, généralement compris entre 70 € pour l’entrée de gamme et 160 ou 170 € pour le haut de gamme, est-il justifié?

Pour ce faire, nous allons comparer les modèles les plus répandus sur le marché Français.

Description du matériel

ModèlePrix (TTC)MicroprocesseurDMIPSFPURAMFlashÉcranPériphériques additionnels
Raspberry pi zero~ 10 €BCM2835 1 Ghz
ARM11 32 bits
1130Oui512 Mo LPDDR2microSD requis
~ 6-7€ pour 16 Go
NonemicroSD
HDMI
GPIO
USB
Wiko Y51 60 €SC7731E
4 Cortex-A7 1.3 Ghz
9880 (multi-core)Oui1 Go8 Go960×4802×5MP camera
Wifi 802 b/g/n
Bluetooth
GSM 2G/3G+
accéléromètre
capteur de lumière
écran capacitif
Carte son
GPS
USB
microSD
Radio FM
TI 82 Advanced~ 70 €Z80 15Mhz
8 bits
0.61Non48 Ko SRAM1 Mo96×64
monochrome
Clavier
USB
Ti 83 premium CE~ 80 €eZ80 48 Mhz
8 bits
7.8Non256 Ko SRAM4 Mo320×240Clavier
USB
Ti Nspire CX II-T CAS~ 150€ARM926EJ-S
396 Mhz
436Oui 64 Mo100 Mo320×240Clavier USB
Casio Graph 35+E II
(fx-9860GIII)
~ 90€SH7305
59 Mhz
106.2Non?? Mo
61 Ko utilisateur
8 Mo128×64
monochrome
Clavier
USB
Casio Graph 90+E
(fx-CG 50)
~ 90 €SH-4A SH7305
117.96 Mhz
212.3Non8 Mo
61 Ko utilisateur
32 Mo396 × 224Clavier
USB
HP Prime G2~165€Cortex A7
528 Mhz
1003Oui256 Mo512 Mo320×240Clavier
USB
NumWorks~ 80 €STM32F730V8T6
Cortex M7
216 Mhz
462Oui256 Ko SRAM8 Mo320×240Clavier
USB

Le tableau ci-dessus décrit divers types de calculatrices ainsi que deux autres matériels: le Raspberry pi zero, un micro-ordinateur vendu sans périphérique (ni écran, ni clavier) et le Wiko Y51, un smartphone Android d’entrée de gamme. Les DMIPS représentent une estimation des performances au benchmark Dhrystone, reflétant les capacités de calcul du microprocesseur. Ces DMIPS ont été critiqués pour être des micro-benchmarks, basés sur une petite quantité de données et de code. Cela va tendre à sous-estimer l’écart entre les processeurs les plus puissant et les moins puissants. Ainsi, le processeur du Raspberry pi zero, avec 1130 DMIPS n’est pas 1130/7.8 = 145 fois plus rapide que la Ti 83 Premium CE. L’écart est certainement beaucoup plus important pour des programmes non triviaux. En réalité, la plupart des applications que l’on fera fonctioner sur un Raspberry pi zero ne pourraient jamais être exécutées sur une Ti 83 Premium CE car elles dépasseraient complètement les capacités mémoire maximales du microprocesseur. Il est aussi à noter que les DMIPS ne prennent pas en compte le calcul des nombres à virgule flottante qui sera extrêmement lent sur la Ti 83 Premium CE car entièrement émulé.

Les caractéristiques techniques du Wiko Y51 dépassent très largement tous les modèles de calculatrices les plus onéreux sur absolument tous les aspects techniques : puissance du microprocesseur, mémoire RAM, mémoire Flash, écran et connectivité. En comparaison à une Ti 82 Advanced, elle a 21845 fois plus de mémoire RAM, un écran qui comporte 75 fois plus de pixels, 8192 fois plus de Flash et un microprocesseur dont les performances sont incommensurables. Son prix est pourtant inférieur.

Du fait de l’absence d’unité de calcul à virgule flottante (FPU), la Ti 83 premium CE doit émuler les nombres à virgule avec les additions et soustractions 8 bits, et utilise une représentation numérique particulièrement peu précise, avec une mantisse de 22 bits, inférieure aux 24 bits d’un FP32 IEEE-754. En bref, cette calculatrice n’est même pas bonne à faire des calculs numériques.

Qu’est-ce qui pourrait justifier ce prix démesuré ?

On ne demande pas forcément d’une calculatrice d’égaler les performances d’un smartphone, mais pourquoi le prix n’est-il pas en rapport avec le matériel ?

Il ne s’agit pas du volume de ventes. Le marché des calculatrices graphiques est très grand, en raison de leur utilité aux examens et concours. C’est un marché international, avec un faible nombre de modèles, puisque Casio, Texas Instruments et Helwet Packard ne vendent chacun que quelques modèles différents. Helwet Packard ne vend plus qu’un seul modèle (HP Prime) alors que Casio vend presque le même matériel à différents prix afin de créer artificiellement une gamme. En comparaison, le Wiko Y51 doit tailler sa place dans un marché encore plus large mais très fragmenté, parmi des milliers de modèles Android différents.

Pourrait-il s’agir du coût du logiciel ? Cela est difficile à dire car les coûts de développement ne sont pas communiqués. Helwet Packard a probablement eu des frais minimes puisqu’il a repris du logiciel libre : FreeRTOS pour le système d’exploitation et GiCalc/Xcas pour le logiciel de calcul formel. Les quelques applications supplémentaires (Classeur, Statistiques) sont triviales à développer. Il est possible que le coût de développement de Texas Instruments soit bien plus élevé, car le microprocesseur Z80 est extrêmement ancien, sous-performant et oblige probablement les ingénieurs à programmer en assembleur en raison des limites de mémoire de la Ti 82 et la Ti 83. Étant donné que le processeur et le logiciel ont faiblement évolué depuis la Ti 81 commercialisée en 1990, il y a quand même trente ans d’amortissement. On remarquera que ce microprocesseur existe depuis 1976 et a représenté, avec le 6502, l’un des deux microprocesseurs 8 bits les plus vendus des années 1980. Le logiciel de Casio évolue aussi lentement, toujours basé sur le Casio Basic qui ne gère que 28 variables nommées de A à Z, plus rho et theta.

Le travail principal des constructeurs semble être la veille réglementaire et des programmes des examens de l’enseignement secondaire des différents pays. Il existe aussi un travail de démarchage auprès des enseignants, qui ont toujours le dernier modèle avant les autres.

Au total, il est probable que les marges bénéficiaires soient très importantes. Cela est d’autant plus problématique qu’il s’agit d’un marché forcé, avec des lycéens, ou plutôt leurs parents, de toute catégorie sociale devant acheter un instrument qui ne servira probablement plus une fois les examens passés.

On peut aussi craindre que cela amplifie les inégalités sociales, en raison d’une gamme de modèles, plus ou moins perfectionnés. Les plus riches pourraient acheter les modèles avec les fonctions les plus avancées, qui leur donnerait un avantage. Ce phénomène pourrait avoir été aggravé par le mode examen, supprimant les options de programmation permettant de débrider les modèles bas de gamme pour leur offrir des fonctionnalités équivalentes aux modèles haut de gamme.

Visite index et cohortes rétrospectives

Problématique

La méthodologie des cohortes rétrospectives n’est presque jamais détaillée dans les articles et pourtant elle est essentielle. Si on s’intéresse, par exemple, au pronostic d’une maladie après diagnostic, on voudra inclure les sujets au moment du diagnostic. Si la recherche s’effectue dans des dossiers médicaux électroniques, les patients atteints de la maladie seront identifiés par une donnée structurée ou non structurée, tel qu’un code diagnostic, le résultat d’un examen biologique ou des mots-clés dans le texte. Il est aussi possible de se baser sur une combinaison de plusieurs éléments (pe. code diagnostic OU biologie). Cette étape conduira à la pré-sélection d’un ensemble de dossiers qui seront ensuite relus, afin de confirmer le diagnostic et les critères d’inclusion.

Trois dates peuvent alors être définies pour chaque sujet:

  1. La date de diagnostic de la maladie (ou date de début des symptômes si c’est ça qui nous intéresse)
  2. La date de la première visite présente dans le dossier médical du centre (ou des centres pour les cohortes multicentriques)
  3. La date index, c’est-à-dire, la date du premier document qui a permis à ce patient d’être pré-sélectionné. Si c’est la combinaison de plusieurs documents (pe. dosage biologique + code diagnostic), alors c’est la date à partir de laquelle les documents sont suffisants pour que le patient ait été pré-sélectionné.

Ces trois dates peuvent différer. Un patient transféré d’un autre centre pourra avoir une date de diagnostic antérieur à la date de la première visite du centre de l’étude. Selon la méthode de pré-sélection, la première visite du centre de l’étude pourra ne pas être identifiée par les codes ou mots-clés recherchés alors qu’une visite ultérieure le sera.

Une erreur méthodologique très fréquente consiste à analyser le pronostic à partir de la date de diagnostic, sans tenir en compte du biais de temps immortel (immortal time bias) induit par le délai entre la date de diagnostic et la date index. Ce phénomène est caricatural si le critère de jugement est la survie globale. Tout sujet décédé avant la date index sera exclu car il ne sera pas pré-sélectionné. S’il y a toujours un intervalle de deux ans entre le diagnostic et la date index, alors on n’observera aucun décès les deux premières années suivant le diagnostic puisque tous les sujets décédés auront été exclus !

Le biais persiste, dans une moindre mesure, si on se base sur la date de la première visite plutôt que la date de diagnostic. Le suivi, en réalité, débute à la date index.

Solutions

Comment peut-on alors correctement modéliser la survie avec les modèles de survie non paramétriques ou semi-paramétriques habituels (Kaplan-Meier et Cox) ?

Date de diagnostic comme baseline

La stratégie idéale, si elle est faisable, consiste à prendre la date de diagnostic comme date de début de suivi mais appliquer une troncature à gauche jusqu’à la date index dans le modèle de survie. Cette troncature à gauche est une fonction assez méconnue. Plutôt que de juste considérer que chaque patient est suivi à partir de T0 jusqu’à une date des dernières nouvelles, faisant alors sortir le sujet de la cohorte, on considère qu’il existe une date des premières nouvelles et une date des dernières nouvelles. Le sujet entre dans la cohorte aux premières nouvelles et en sort aux dernières nouvelles. Le nombre de sujets à risque peut alors croître puis décroître, puisqu’il y a des gagnés de vue et des perdus de vue. Cette méthode permet de conserver la forme de la courbe de survie originale, en reposant sur l’hypothèse d’entrée et de sortie de la cohorte au hasard, ou, du moins, pour des raisons non corrélées à l’outcome.

Méthode landmark

Cette stratégie a une limite: elle nécessite qu’un nombre suffisant de sujets ait une date index égale à la date de diagnostic, sinon, dans le pire des cas, on commencera par un estimateur de Kaplan-Meier à 0/0, ce qui rendra impossible toute estimation de courbe de survie. Si la majorité voire la totalité des sujets ont un écart important entre la date de diagnostic et la date index (pe. 1 an), on peut comprendre que l’évolution initiale est impossible à connaître. On doit alors renoncer à la comparaison sur la période initiale de la courbe de survie. On peut utiliser la méthode landmark, qui consiste à redéfinir la baseline, c’est-à-dire, le T0 de la courbe de survie, au diagnostic+constante tel qu’un an après le diagnostic. Au nouveau point de départ, tous les sujets auront la même ancienneté de la maladie. Ils seront suffisamment nombreux pour que le tracé de la courbe soit possible. Tout sujet ayant eu l’événement avant le landmark sera exclu (censure à gauche).

Date index comme baseline

Une stratégie alternative consiste à définir la baseline (T0 de la courbe de survie) comme la date index. Cette méthode fournit la meilleure précision statistique, car garantit un échantillon de taille maximale à T0. Par contre, cette méthode fournit des courbes de survie d’allure exponentielle quand bien même ça ne reflète pas du tout l’évolution de la maladie. En mélangeant tous les stades d’ancienneté de la maladie à baseline, le rythme d’apparition des événements devient une moyenne des risques associés à chaque ancienneté. Par exemple, l’ataxie spinocérébelleuse de type 2 est une maladie neurologique dégénérative d’évolution progressive lente mais inexorable. La figure 1 de l’article « Prediction of Survival With Long-Term Disease Progression in Most Common Spinocerebellar Ataxia » (doi: 10.1002/mds.27739, PMID: 31211461) décrit une survie globale de 97.8% à 10 ans, 78% à 20 ans, 31% à 30 ans et 11.2% à 40 ans. Ainsi, la mortalité précoce (< 10 ans) est négligeable alors l’issue fatale survient majoritairement entre 15 et 35 ans. On observe pourtant des courbes très différentes sur la figure 1 de l’article intitulé « Survival in patients with spinocerebellar ataxia types 1, 2, 3, and 6 (EUROSCA): a longitudinal cohort study » (doi: 10.1016/S1474-4422(18)30042-5, PMID: 29553382). Le taux de survie à 10 ans du l’ataxie spinocérébélleuse de type 2 est environ de 73%, avec un rythme de décès semblant assez constant sur les 10 années de suivi. Cela est explicable par un T0 correspondant à la date index, et concernant des cas prévalents d’ancienneté très variable. Le rythme de décès est alors égal à la moyenne des rythmes de décès de toutes les anciennetés, pondérée par la prévalence des anciennetés. Cette attitude peut néanmoins se défendre pour les situations où l’évolution est peu dépendante de l’ancienneté, c’est-à-dire, correspondant à des courbes de survie d’allure exponentielle.

Méthodes paramétriques ?

Éventuellement, on pourrait aussi utiliser la loi de Weibull pour compléter le début du suivi mal connu.

Généralisation des concepts

La distinction entre date de début des symptômes et date de diagnostic peut parfois avoir une importance (pe. syndrome démentiel) et c’est souvent la date de début des symptômes qui importe le plus, avec néanmoins des problèmes de difficulté à mesurer la date de début des symptômes.

Même s’il est fait référence au dossier médical informatisé, les problèmes décrits dans ce billet sont tout aussi applicables aux dossiers papiers, pour lesquels il existe des documents ou source de données index. Il est aussi généralisable aux registres. Il est enfin généralisable aux cohortes prospectives incluant des cas prévalents, pour lesquels la visite index sera généralement la visite d’inclusion.

Concept apparenté, l’anti-cohorte ou cohorte inversée

Ne cherchez pas ce concept dans la littérature scientifique, vous ne l’y trouverez pas. Le terme de cohorte inversée ou anti-cohorte est une invention de l’auteur de ce blog. Il s’agit d’une méthodologie à classer dans la sémiologie fongique des études, c’est-à-dire, un exemple de ce qu’il ne faut pas faire. Plutôt que de sélectionner les sujets nouvellement diagnostiqués sur une période donnée (pe. entre 2010 et 2019) et de les suivre jusqu’à survenue d’un événement ou censure administrative à une date de point (pe. 31 décembre 2019), la cohorte inversée sélectionne les sujets sur la présence d’un événement sur une période récente (pe. 2018-2019), puis on remonte le dossier médical jusqu’à la visite la plus ancienne du dossier médical (pe. 2000-2019) et on l’analyse comme si le sujet avait été sélectionné sur cette première visite.

Le taux de survenue d’événement atteint alors 100%. Cela ressemble à une méthodologie cas-témoin dans laquelle il n’y aurait que des cas. Malheureusement, le fichier de données se présente comme celui d’une cohorte et un statisticien non au fait de la méthodologie de sélection employée pourrait alors l’analyser comme s’il s’agissait d’une cohorte ordinaire. Les sujets dont la date de première visite est récente auront alors une survie avant événement raccourcie.

Une variante de cette cohorte inversée est la sélection de la « file active » des patients (pe. ayant eu une visite entre janvier et décembre 2019) qu’ils aient eu ou non l’événement d’intérêt, puis de remonter jusqu’à la première visite du dossier médical (pe. jusqu’à l’an 2000 pour certains patients). La visite index est alors en 2019, avec une forte représentation des cas prévalents. Même si certains événements peuvent survenir après la date index (pe. date index en février 2019 et événement en septembre 2019), si on utilise la date de première visite comme baseline du suivi, on aura une période d’immortalité extrêmement longue et un très faible nombre d’événements au total. Toute variable positivement corrélée à la date de première visite sera corrélée négativement corrélée au délai avant survie.