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:

La comparaison brute des scores post-intervention, sans ajustement (p.e. test de Student sur séries indépendantes)
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.
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:

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
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 :

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.
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.