Conditions de validité du coefficient de corrélation de Pearson et de la régression linéaire

Pour faire court, la seule condition de validité pour le calcul d’un coefficient de corrélation de Pearson ou l’estimation d’une régression linéaire, est l’existence d’une variance non-nulle sur chacune des deux variables, sous peine de division par zéro. Tous les autres éléments ne sont que des « précautions d’interprétation ».

La régression linéaire, estimée par les moindres carrés a des propriétés générales très intéressantes, parce l’estimateur des moindres carrés, c’est celui de la moyenne. Cela garantit une espérance d’erreur nulle, là où la plupart des estimateurs de modèle n’ont pas cette propriété extrêmement utile, voire indispensable dans certaines situations.

Le coefficient de corrélation de Pearson peut avoir une interprétation conjointe à la régression linéaire simple (son carré étant la variance expliquée par la régression) mais bénéficie aussi de propriétés propres très intéressantes.

Ce billet lance une réflexion sur l’usage de ces statistiques afin d’aller au-delà des algorithmes automatiques tournées autour de l’hypothèse de normalité ; cette dernière étant aussi plausible que l’existence du monstre du Loch Ness comme le montre cet algorithme.

Définition du coefficient de corrélation de Pearson

Le coefficient de corrélation de Pearson entre deux variables est égal à l’espérance du produit des variables préalablement centrées-réduites. C’est aussi égal à la covariance de ces deux variables centrées-réduites.

Cette définition, n’implique absolument pas d’hypothèse de normalité, de linéarité, d’homoscédasticité ou quoi que ce soit du genre. C’est juste un calcul mathématique, qui lui confère quelques propriétés générales très précieuses.

Propriétés générales du coefficient de corrélation de Pearson

Considérant deux variables A et B, de lois quelconques, dépendantes ou pas, alors

COV(A,B)=COR(A,B)×sqrt(VAR(A)×VAR(B))

Où VAR représente la variance, COV la covariance et COR, la corrélation de Pearson.

Ainsi, connaissant, l’écart-type (ou la variance) de A et de B, on peut calculer la variance de la somme de A et B ou de la différence, en s’aidant des deux formules

VAR(A+B) = VAR(A)+VAR(B)+2×COV(A,B)

VAR(A-B) = VAR(A)+VAR(B)-2×COV(A,B)

Ces propriétés mathématiques ne sont pas asymptotiques, elles sont exactes, même sur un échantillon fini. Même avec n=3 observations. Il n’y a aucune approximation et aucune hypothèse distributionnelle.

C’est une conséquence directe de la décomposition de la formule (a+b)²=a²+b²+2×a×b.

Ces propriétés sont très utiles pour les méta-analyses et/ou l’extraction de données d’un article.

Considérons la situation assez fréquente dans laquelle on souhaite connaître la variance intra-sujet d’une mesure afin de calculer le nombre de sujets nécessaires dans une étude en cross-over. À partir d’un article fournissant la moyenne et l’écart-type de la mesure d’intérêt, dans la population d’intérêt, à deux temps différents, et fournissant le petit p d’un Student sur séries appariées, on peut remonter du petit p à la statistique t de Student, pour en déduire la variance de la différence. On peut aussi aller plus loin, en utilisant la variance de la différence et la variance à chaque temps pour en déduire la covariance entre les deux temps ainsi que le coefficient de corrélation de Pearson. Cela fonctionne sans hypothèse de normalité. Quand bien même le petit p du Student sur séries appariées a une validité dépendant d’une hypothèse de normalité asymptotique, le calcul réalisé pour remonter au coefficient de corrélation reste vrai même si l’hypothèse est fortement fausse.

Si les auteurs omettent le test de Student sur séries appariées mais fournissent un coefficient de corrélation de Pearson, alors grâce aux moyennes et aux écarts-types, on peut calculer la variance d’une différence appariée, ainsi que faire le test de Student sur séries appariées. En bref, variance, covariance, corrélation de Pearson ont des propriétés mathématiques les rendant complémentaires, permettant l’extraction de données et la vérification de cohérence. Soit dit en passant, cela fonctionne tout aussi bien avec des variables binaires qu’avec des variables quantitatives continues ou discrètes.

Sur la population, le coefficient de corrélation de Pearson entre deux variables X et Y est égal au coefficient de corrélation de Pearson entre deux moyennes Mx et My d’échantillons de taille n d’observations indépendantes et identiquement distribuées, et ce pour tout n. Cela se démontre, en passant par la covariance:

COR(Mx, My) = COV(Mx, My)/sqrt(VAR(Mx)×VAR(My)) = (COV(X,Y)/n)/sqrt(VAR(X)/n × VAR(Y)/n) = COV(X,Y)/sqrt(VAR(X)×VAR(Y)) = COR(X,Y)

Au fond, c’est cette propriété qui permet d’estimer une différence de moyennes appariée en connaissant seulement de la moyenne et l’écart-type estimés sur chaque série de données ainsi que d’une estimation du coefficient de corrélation entre les deux variables.

Le code suivant montre qu’en présence d’un effet non-linéaire, le coefficient de corrélation de Pearson vérifie bien cette propriété alors que le coefficient de corrélation de Spearman entre les moyennes est dix fois plus fort (R=0.70) qu’entre les variables X et Y (R=0.07).

set.seed(2020)
a=rexp(1e6)
b=2*(a-1)^2+rexp(1e6)
cor.test(a,b, method="spearman") # Spearman's r=0.07
cor.test(a,b, method="pearson") # Pearson's r=0.70

m=sapply(1:1e6, function(x) {
	a=rexp(100)
	b=2*(a-1)^2+rexp(100)
	c(mean(a), mean(b))
})

cor.test(m[1,], m[2,], method="spearman") # Spearman's r=0.70
cor.test(m[1,], m[2,], method="pearson") # Pearson's r=0.70

Malheureusement, un coefficient de corrélation de Spearman entre les variables X et Y ne permet absolument pas d’estimer la corrélation de Spearman entre les variables Mx et My sauf si on fait des hypothèses distributionnelles fortes (p.e. relation linéaire et résidus normaux), et même ainsi, la relation n’est pas une simple égalité.

Définition théorique de la régression linéaire

Le modèle de régression linéaire simple sur un échantillon de taille n comportant des variables (X1, …, Xn) et (Y1, …, Yn) est généralement noté :

Yi=beta0+beta1×Xi+Epsilon_i

Pour i=1, …, n

où beta0 et beta1 sont des constantes.

On suppose que les Epsilon_i suivent tous une loi normale, de même variance et que tous les Epsilon_i sont indépendants. On associe généralement à ce modèle, l’estimateur des moindre carrés.

Que faire lorsqu’on sait qu’aucune loi du vivant ne suit une loi normale, comme décrit dans ce billet ? Cela rend-il impossible d’utiliser ce modèle ?

Ma réponse, c’est qu’en redéfinissant ce modèle, il apparaît que ces conditions ne sont pas nécessaires.

Redéfinition de la régression linéaire

Cette section est un peu technique et pas nécessaire à la compréhension globale du billet. N’hésitez pas à passer à la section « propriétés de la régression linéaire » si vous ne comprenez rien. Elle me permet juste de répondre aux puristes qui diraient que l’homoscédasticité, la linéarité et la normalité des résidus font partie, par définition, du modèle linéaire et qu’en conséquence tout ce que je dis n’est qu’un tas d’ineptie puisque je remets en cause une définition. Je suis peut-être hérétique en remettant en cause un dogme une définition.

L’estimateur des moindres carrés de la régression linéaire est applicable à n’importe quel échantillon (sauf division par zéro dans des cas exceptionnels). Il permet d’estimer des coefficients beta0 et beta1 ainsi que des prédictions comme des combinaisons linéaires de beta0 et beta1. Ce sont des statistiques, soumises à des fluctuations d’échantillonnage.

À partir d’une statistique calculable sur un échantillon, on peut toujours définir une statistique dans une population infinie, comme la limite, quand n tend vers l’infini, de l’espérance de la statistique sur un échantillon aléatoire de taille n, d’observations indépendantes et identiquement distribuées issues de la population considérée. Il est théoriquement possible que cette limite n’existe pas, comme avec la loi de Cauchy, mais dans les cas usuels en biologie et médecine où les distributions sont bornées, et avec la régression linéaire qui a des bonnes propriétés de stabilisation asymptotique, elle devrait toujours exister.

Ainsi, en transformant l’estimateur en statistique dans la population, je définis les coefficients (beta0 et beta1) et prédictions d’un modèle linéaire de manière totalement indépendante de la moindre hypothèse de linéarité, homoscédasticité ou normalité. La seule hypothèse est l’existence d’une limite à ces coefficients. Cela me permet alors d’analyser les propriétés de ces statistiques dans le cas ordinaire de non-respect des hypothèses théoriques de validité.

Propriétés de la régression linéaire

La moyenne est une statistique directement liée à l’estimateur des moindres carrés. En effet, c’est la statistique qui minimise les carrés des écarts à cette statistique.

Vous pouvez le vérifier par vous-même sous le logiciel R:

v=c(1,30,50)
mean(v) # vaut 27
optimize(function(position) {sum((v - position)^2)}, interval=range(v)) # vaut 27
coef(lm(v ~ 1)) # vaut 27

La régression linéaire est aussi habituellement estimée par les moindres carrés. Cela lui confère des propriétés intéressantes, sans hypothèse de normalité, homoscédasticité ni linéarité.

Sur un échantillon fini, aussi petit soit-il, la moyenne des résidus d’une régression linéaire estimée par les moindres carrés est nulle. Les résidus positifs sont compensés par les résidus négatifs. Cela est toujours vrai, sans la moindre hypothèse (même pas l’indépendance des observations). Le code R suivant permet de le vérifier:

x=rexp(3)
y=x^(1/3)+rexp(3)
mean(resid(lm(y ~ x))) # zéro, aux erreurs d'arrondi des nombres à virgule flottante près

Si on considère une régression linéaire des moindres carrés dont les coefficients sont connus exactement (ou estimés sur un très grand échantillon représentatif de la population), l’espérance des erreurs de prédiction de ce modèle est nulle sur cette même population, sous l’hypothèse d’observations indépendantes et identiquement distribuées. Les erreurs positives compensent exactement les erreurs négatives. Cela est faux avec d’autres modèles linéaires tel que le modèle linéaire identité-binomial estimé par le maximum de vraisemblance. Vous pouvez vérifier cette propriété avec le code R suivant:

set.seed(2020)
x=rexp(1e7) # distribution non normale
population=data.frame(x=x, y=x^(1/3)+rexp(1e7)) # relation non linéaire
model = lm(data=population, y ~ x) # coefficient presque exactement connus
x=rexp(1e6)
bigsample=data.frame(x=x, y=x^(1/3)+rexp(1e6))
mean(predict(model, newdata=bigsample) - bigsample$y) # erreur de prédiction moyenne presque nulle

Prenons l’exemple d’une mutuelle santé à but non lucratif voulant calculer les coûts de cotisation en adaptant le coût de cotisation à l’âge, qui est une variable reflétant la consommation de soins. Pour chaque cotisant, on peut calculer un coût de consommation précis (en euros). Les frais de gestion sont éventuellement ventilés proportionnellement au coût de cotisation, de telle sorte qu’on puisse calculer un « coût » de chaque cotisant de telle sorte que les dépenses totales de la structure sont égale à la somme de tous les coûts de tous les cotisants. L’échantillon exhaustif des cotisants sur les trois années précédentes sert ensuite à estimer les paramètres d’une régression linéaire des moindres carrés expliquant le coût par l’âge, permettant de calculer les cotisations mensuelles, adaptées à l’âge, pour l’année suivante. À moins que la pyramide des âges ou le profil de consommation selon l’âge évolue beaucoup l’année suivante, cette méthode garantit l’équilibre budgétaire, quand bien même la relation entre la consommation et l’âge n’est pas linéaire. On peut noter que la mutuelle peut aussi préférer un système plus solidaire, avec une cotisation indépendante de l’âge, en supprimant l’âge du modèle, conduisant alors à une cotisation égale à la moyenne des coûts des cotisants. Cette moyenne conserve aussi la propriété d’équilibre budgétaire. L’usage de modèles « robustes », excluant par exemple, le pourcent de consommateurs les plus forts, ne garantirait pas l’équilibre budgétaire.

Même si l’exemple fourni est basé sur une régression linéaire simple, cette bonne propriété d’espérance d’erreur nulle est applicable au modèle linéaire général, et donc à des modèles multivariés prenant en compte plusieurs paramètres pour le calcul des cotisations. Même en présence d’interactions, cette propriété est conservée.

Relation entre régression linéaire et coefficient de corrélation de Pearson

La régression linéaire des moindres carrés minimise la variance résiduelle. La variance expliquée, ou R², est égale à un moins le rapport entre la variance résiduelle et la variance totale. Pour une régression linéaire simple, le coefficient de corrélation de Pearson élevé au carré est égal à ce R². Cela est toujours vrai sur un échantillon, encore une fois, sans hypothèse particulière.

Propriétés spécifiques du modèle linéaire

Hypothèse de linéarité respectée

Définissons l’existence d’une relation linéaire entre une variable Y et une variable X par l’existence d’un modèle linéaire tel que l’espérance de la variable Y conditionnelle à une valeur de X est égale à la prédiction de Y par X. L’estimateur des moindres carrés est capable de trouver les bons coefficients, garantissant cette propriété, quand bien même il y a une hétéroscédasticité et les distributions de Y conditionnelles à chaque valeur de X diffèrent en forme les unes des autres ! On repose seulement sur l’hypothèse d’observations indépendantes et identiquement distribuées dans la distribution bivariée (X,Y). Le code R suivant illustre cette propriété:

set.seed(2020)
x=rep(c(1,2,3), c(1e5, 2e6, 1e6))
y=c(2+rnorm(1e5), # résidu normal de faible variance
4+(rexp(2e6)-1)*3, # résidu selon loi exponentielle d'écart-type égal à 3
6+(runif(1e6)-0.5)*10 # résidu selon loi uniforme, de très forte variance
)

predict(lm(y ~ x), newdata=data.frame(x=1:3)) # prédit bien 2, 4, 6

Ainsi, avec des coefficients estimés par l’estimateur des moindres carrés sur un grand échantillon, l’erreur de prédiction de Y conditionnelle à n’importe quelle valeur de X, est nulle. Cela garantit notamment que l’espérance de l’erreur de prédiction reste nulle quand bien même la distribution de X change.

On remarquera que l’estimateur des moindres carrés avec une quelconque pondération aura la même espérance des coefficients. Par exemple, une pondération ayant pour objectif de « rectifier » les problèmes d’hétéroscédasticité, sera asymptotiquement équivalent à la version non pondérée.

Hypothèse d’indépendance entre les résidus et la variable X respectée

De manière équivalente, on peut dire que la distribution de la différence entre la prédiction et l’observation est indépendante de X, c’est-à-dire, la distribution de cette différence conditionnelle à une valeur de X a la même distribution quelle que soit la valeur de X.

Si cette propriété est vérifiée, alors, non seulement on garantit que conditionnellement à chaque valeur de X, l’espérance de l’erreur est nulle, mais on peut même prédire la distribution de cette erreur. On peut empiriquement déterminer cette distribution comme la distribution des résidus observés. Il est possible de lisser, ou non, la distribution empirique de ces résidus par une estimation par noyau. En s’aidant de la distribution empirique de X estimée par noyau, on peut fournir une estimation de la distribution jointe (X,Y) avec densité de probabilité bivariée. Néanmoins, je déconseille fortement ce type de modélisation car cette hypothèse est généralement (toujours?) fausse. On peut souvent utiliser directement la distribution bivariée (X,Y) empirique, plus ou moins lissée si nécessaire. C’est néanmoins plus intéressant lorsqu’on s’intéresse à la distribution de Y conditionnelle à une valeur précise de X puisque dans ce cas, la distribution de Y conditionnelle à X est susceptible de n’être estimable qu’avec zéro ou une valeur si on souhaite une estimation empirique sans hypothèse de modélisation.

Hypothèse de normalité des résidus

Si cette hypothèse est vérifiée (ce qui n’arrive jamais ?), alors… elle est vérifiée et on peut reposer dessus. Cela veut dire, par exemple, qu’on peut estimer la variance résiduelle empiriquement, puis utiliser cette variance résiduelle comme paramètre d’une distribution normale afin de connaître la distribution de l’erreur résiduelle conditionnelle à n’importe quelle valeur de X. Par rapport à l’hypothèse précédente, on peut gagner un peu en précision sur l’estimation de cette distribution conditionnelle sur des échantillons de taille modeste. Reposer sur cette propriété engendre un biais mais est susceptible de réduire l’erreur sur de petits échantillons, lorsque l’écart à la normalité est modeste, par rapport à l’estimation empirique de la distribution des résidus. C’est alors un choix guidé par le rapport biais/erreur.

Normalité asymptotique

Si vous avez déjà essayé d’estimer des paramètres de régression Passing BaBlok par bootstrap non paramétrique sur un échantillon de taille modeste, vous avez dû remarquer que les fluctuations d’échantillonnages du Passing BaBlok sont discrètes. C’est dû au fait que la distribution empirique est discrète sur un échantillon de taille modeste. Le Passing BaBlok a des fluctuations d’échantillonnages chaotiques lorsque X ou Y suivent des lois discrètes. Ce problème n’existe pas avec la régression linéaire des moindres carrés dont les estimateurs de coefficients suivent asymptotiquement une loi normale multivariée quelque soient les distributions de X et Y à condition que l’échantillon soit constitué d’observations indépendantes identiquement distribuées. Cette normalité asymptotique s’applique aussi au coefficient de corrélation de Pearson. Il n’y a pas à reposer sur une hypothèse d’homoscédasticité ou de normalité des résidus.

Intervalles de confiance et petits p

Sur des échantillons de taille suffisante, alors, le boostrap permet de fournir des intervalles de confiance asymptotiquement corrects aussi bien pour la régression linéaire que pour le coefficient de corrélation de Pearson. Les estimateurs d’intervalles de confiance classiques peuvent par contre être biaisés. La transformation z de Fisher et l’approximation à une loi normale du coefficient de corrélation de Pearson sont asymptotiquement corrects mais l’estimateur de sa variance comme égal à 1/(n-3) est parfois asymptotiquement biaisé. Cette approximation est asymptotiquement correcte lorsque le vrai coefficient de corrélation nul, mais d’une manière générale elle est susceptible d’être asymptotiquement biaisée. De même, l’intervalle de confiance de Wald sur les coefficients d’une régression linéaire est susceptible d’être asymptotiquement biaisé. Si on a des résidus normalement distribués, indépendants et identiquement distribués, alors les approximations sont correctes. C’est pourquoi je recommande l’usage du boostrap dans le cas général.

Sensibilité aux outliers

L’estimateur des moindres carrés et fortement influencé par les valeurs atypiques (outliers) et donc la régression linéaire comme le coefficient de corrélation de Pearson le sont aussi. Selon les cas, cela peut-être souhaitable ou pas. Reprenons l’exemple d’une mutuelle santé qui s’intéresse au coûts de ses prestations. Si on s’intéresse aux rentrées d’argent associés à chaque adhérent, exprimées comme la somme des cotisations moins les dépenses associées aux remboursements, alors la majorité des adhérents fourniront une rentrée d’argent positive alors que quelques rares adhérents coûteront des dizaines de milliers d’euros en frais d’hospitalisation. Un adhérent qui coûte 70 000 €, pèse 700 fois plus sur le budget, qu’un adhérent qui coûte 100 €, et cela doit forcément être pris en compte si on souhaite l’équilibre budgétaire. En conséquence, si on veut garantir l’équilibre budgétaire, il faut une estimation précise de la proprotion de sujets qui coûtent 70 000 €. Vous comprenez bien que ce n’est pas avec 30 observations qu’on pourra estimer une moyenne correcte ! La règle selon laquelle la méthode de Student fonctionne dès que n>= 30 est ridicule puisque ça dépend fondamentalement de la fréquence et le degré d’atypie des valeurs atypiques. Ce phénomène est encore plus marqué pour les assurances qui remboursent des frais en cas d’accident très onéreux et très rare. La distribution des risques doit être alors estimés par des modèles bien plus complexes qu’un échantillon avec n=30. Il est bien évident que l’espérance reste le paramètre clé. Il ne faut surtout exclure ces outliers ou se baser sur la médiane. La sensibilité aux outliers est une nécessité; ce sont eux qui contiennent l’information.

Les choses sont très différentes si on s’intéresse, par exemple, à la corrélation entre la vitesse de sédimentation (VS) et la C Reactive Protein (CRP). L’espérance n’a plus d’importance et on s’intéressera plutôt à l’idée de seuils de positivité ou à des seuils pathologiques. Il paraîtra toujours pertinent de considérer qu’une valeur de CRP à 300 est supérieure à une valeur à 50, mais le ratio 300/50=6 ne revêtira pas de pertinence en pathogenèse. Devant cette situation, la corrélation sera mieux appréciée par un tau de Kendall ou un coefficient de corrélation de Spearman que par un coefficient de corrélation de Pearson. La régression linéaire posera des problèmes de stabilité d’estimateur sur des échantillons de taille modeste, mais manquera aussi de pertinence dans la description des relations. Le fait de perdre en performance prédictive sur les valeurs typiques pour améliorer la description des outliers, pourrait être contre-productif. Cela dépend néanmoins de l’usage de cette relation. Un meilleur exemple pourrait être la description de la relation entre deux techniques de dosage des anticorps ciblant un même antigène. Une modèle linéaire ou non-linéaire avec un estimateur robuste aux outliers pourrait être utilisé pour convertir un dosage en l’autre et établir ainsi une équivalence.

Résistance aux distributions discrètes

Juste en passant, le coefficient de corrélation de Spearman n’est pas spécialement adapté aux distributions discrètes. Le coefficient de corrélation de Pearson est parfaitement calculable avec des distributions discrètes, voire binaires, et ne souffre pas d’instabilité parce qu’il n’y a généralement pas d’outliers dans ce contexte.

Au contraire, le coefficient de corrélation de Spearman, avec sa transformation des rangs, va créer des écarts entre deux valeurs successives, d’autant plus grandes que la valeur est fréquente, rendant plus délicate l’interprétation du coefficient, sans compter les problèmes d’estimation de sa variance. On peut toujours s’en tirer avec du bootstrap si l’échantillon est de taille suffisante.

La régression linéaire est tout à fait pertinente sur des variables binaires. Si Y et X sont toutes deux binaires, la pente de la régression linéaire s’interprète comme la différence de proportions de Y=1 entre le groupe où X=0 et le groupe où X=1. L’ordonnée à l’origine (intercept) s’interprète comme la proportion de Y=1 dans le groupe où X=0.

Limites d’interprétation

En cas de relation non linéaire, il peut exister une corrélation très forte entre deux variables, mais le coefficient de corrélation de Pearson peut-être nul (ou très faible) et la régression linéaire peut avoir une pente nulle (ou très faible) et une très faible variance expliquée:

a=rnorm(1e4)
b=a^2
cor(a,b) # corrélation de Pearson nulle
coef(lm(b~a)) # pente nulle
plot(a,b)

Dans ces conditions, des modèles non linéaires permettent de prédire la valeur d’une variable en fonction de l’autre alors que le modèle linéaire n’a pas plus de pertinence que de fournir une moyenne générale. Le coefficient de corrélation de Pearson garde toutes ses propriétés intéressantes, mais ne peut pas être interprété comme une force d’association entre les deux variables. Il peut toujours s’interpréter comme la racine carrée de la variance expliquée par le modèle linéaire, qui est alors presque nul.

En cas d’hétéroscédasticité, le modèle linéaire reste toujours pertinent mais on peut espérer une meilleure stabilité des estimations en pondérant les résidus afin de rectifier l’homoscédasticité. Le modèle a toujours une erreur de prédiction moyenne nulle, voire une erreur de prédiction moyenne nulle conditionnellement à toute valeur de X, si la relation est linéaire, mais la distribution exacte des résidus diffère selon la valeur de X et ne peut pas être juste calculée comme la distribution empirique des résidus.

En cas de non linéarité, le modèle linéaire garde toujours la propriété d’erreur moyenne nulle, mais la variance résiduelle est susceptible d’être bien plus élevée que dans un modèle linéaire, de telle sorte qu’on peut dire qu’il a une performance prédictive médiocre. On perd aussi la propriété d’espérance nulle de l’erreur conditionnelle à toute valeur de X. Enfin, un changement de distribution de X, n’affectant pourtant pas la relation, peut faire apparaître un biais d’estimation, c’est-à-dire, une espérance d’erreur non nulle.

Conclusion

Beaucoup de modèles reposent sur une structure et des hypothèses manifestement fausses. On n’utilise jamais vraiment ces modèles (heureusement car ils sont faux) mais seulement leurs estimateurs. Il me paraît important d’étudier le bon ou mauvais comportement de ces estimateurs dans les cas usuels, en ne faisant qu’un minimum d’hypothèses. J’ai montré ainsi que l’estimateur des moindres carrés du modèle linéaire a des propriétés générales très intéressantes.

Le coefficient de corrélation de Pearson n’est pas un modèle mais une statistique, ayant de nombreuses propriétés intéressantes, comme j’ai montré plus haut. Il n’a pas vraiment de condition de validité mais seulement des limites d’interprétation dans certaines situations.

Le petit p bidon sous-puissant

Quelques définitions sur les risques

Considérant que la plupart des études font des comparaisons bilatérales, qu’il s’agisse d’épidémiologie ou de recherche clinique, l’hypothèse nulle est généralement l’absence totale d’effet de l’intervention ou l’exposition considérée. La plausibilité de cette hypothèse nulle est généralement douteuse, notamment pour les essais cliniques dans lesquels la question porte plus sur le signe et l’amplitude de l’effet que sur la réalité d’un effet. En bref, il paraît absolument impensable que le traitement médical ou chirurgical de la hernie discale ait exactement le même résultat fonctionnel à 1 an. Quand je dis, exactement le même résultat fonctionnel, c’est que la moyenne d’une échelle fonctionnelle serait identique même avec 10 000 chiffres après la virgule. Les vraies questions qui se posent sont :

  1. Lequel des deux traitements est le meilleur (signe de la différence) ?
  2. Est-ce que les deux traitements ont un résultat fonctionnel moyen presque équivalent ou très différent ?

Au sens strict, c’est la balance bénéfices/risques qu’on doit évaluer, en prenant en compte les effets indésirables médicamenteux et les complications chirurgicales, mais pour simplifier, on se concentre sur le résultat fonctionnel dans un premier temps.

Considérant donc que les deux traitements ne peuvent pas être strictement équivalents, le risque l’erreur de première espèce n’existe pas, au sens strict, mais elle est remplacée par les deux risques suivants :

  1. Conclure à l’existence d’une différence dans le sens opposé de la réalité (erreur de troisième espèce).
  2. Conclure à l’existence d’une différence importante alors que la différence réelle est totalement négligeable (quasi-équivalence). J’appellerai ça l’erreur de type Ib.

Définition du petit p bidon sous-puissant

Lorsqu’on utilise un échantillon trop petit et qu’on recherche un effet modeste, alors le signal (effet réel) devient négligeable par rapport au bruit (erreur aléatoire), de telle sorte qu’avec un seuil de significativité bilatéral à 5%, on a 2.5% de chances de conclure à une différence significative dans un sens et 2.5% de chance de conclure à une différence significative dans l’autre sens.

Dans ce contexte, le petit p est indépendant de la différence réelle. Un petit p significatif n’apporte plus aucune information. Il n’aide pas à identifier la réalité d’une différence non négligeable puisqu’il a la même probabilité d’arriver que la différence soit nulle, négligeable, ou non négligeable. Il n’aide pas non plus à identifier la direction de la différence puisqu’il a autant de chances d’aller dans le bon sens que dans le sens opposé et est donc indépendant du signe de la différence et donc non informatif dessus.

Je parlerai donc de petit p bidon sous-puissant pour décrire les petits p significatifs dans une situation de rapport signal/bruit très proche de zéro. Cela regroupe donc, trois cas :

  1. La différence réelle est négligeable ou nulle, mais le petit p est significatif
  2. La différence réelle n’est pas négligeable mais le petit p va dans le mauvaise direction (erreur de troisième espèce)
  3. La différence réelle n’est pas négligeable et le petit p va dans la bonne direction, mais il y avait autant de chances que ça aille dans la direction opposée, de telle sorte, que ce n’est que par pure chance que la conclusion de l’étude est correcte.

Considérer que le troisième item est bidon peut vous choquer, mais cela me paraît pertinent du point de vue de l’information. De mon point de vue, des propos peuvent être considérés de bidon, s’ils sont totalement indépendants de la réalité, ce qui implique de dire parfois des choses vraies et parfois des choses fausses. Quelque chose est bidon à partir du moment où il est décorrélé de la réalité. À l’opposé, dire systématiquement le contraire de la réalité, c’est informatif, puisqu’on peut alors croire le contraire de ce qui est dit.

Sémiologie

Devant un petit p significatif, certains signes évoquent un petit p bidon sous-puissant.

  1. Échantillon de toute petite taille face à l’effet attendu, suggérant une puissance très faible
  2. Estimation ponctuelle démesurée en comparaison à ce qui paraît plausible
  3. Petit p à la limite de la significativité (p typiquement compris entre 0.01 et 0.05)
  4. Multiplicité des tests, apparente ou cachée
  5. Autres tests répondant à la même question ne montrant pas plus de tendance à la significativité que ce qui est explicable par le hasard (environ un petit p sur 20 significatif, la moitié du temps dans le sens opposé à ce que veulent montrer les auteurs)
  6. Lorsqu’une différence d’évolution (p.e. Student inter-groupe sur différences intra-sujets) est analysée sur un paramètre dont la stabilité est attendue, apparition d’une dégradation importante dans le groupe contrôle et d’une amélioration de même amplitude dans le groupe expérimental et une différence à baseline qui va dans le sens opposé à la différence finale.
  7. Lorsqu’une différence d’évolution est analysée sur un paramètre qui doit évoluer, la moitié de la différence des différences est due à la différence à baseline et l’autre moitié à la différence finale.

Certains de critères sont très subjectifs, telle que les deux premiers, mais lorsque beaucoup d’éléments sont présent et fortement marqués, le tableau est évocateur.

À un niveau plus global de la littérature, une méta-analyse peut estimer de manière à peu près correcte l’effet réel avec des milliers de patients, ce qui permet ensuite de mieux identifier les études sous-puissantes. On peut craindre un biais de publication ainsi qu’un selective reporting biais dans ces études.

La multiplicité des tests cachée a sa propre sémiologie:

  1. Écart au protocole, sur les analyses (suspect lorsque celui-ci est disponible)
  2. Critères de jugements présentés dans la partie méthodes mais pas dans les résultats
  3. Grande majorité (voire totalité) de petits p significatifs (suggérant une très bonne puissance sur tous les tests réalisés) dans les résultats, mais presque tous compris entre 0.01 et 0.05 (alors qu’en cas de puissance à 90%, on en a 50% de petits p en dessous de 0.0012)
  4. Étrangeté des analyses qui « tirent dans les coins », comme la corrélation entre le max d’un dosage biologique entre J1 et J4 corrélé à la mortalité en réa, puis la moyenne d’un autre dosage biologique entre J2 et J3 corrélé à la mortalité intra-hospitalière, tous deux dans le même article.
  5. Critère de jugement principal inattendu étant donné la population et l’intervention, voire disparition des critères de jugements attendus.

Explication des critère N°6 et 7

Je crains que la pertinence de ce critère ne paraisse pas être une évidence au novice. Pour le comprendre, il faut raisonner en termes de distributions conditionnelles au petit p significatif. Pour cela, je partirai d’un cas d’école. Nous allons expliquer le critère N°6. Le N°7 est une variante assez simple qui en découle directement.

L’article est intitulé « Special nutrition intervention is required for muscle protective efficacy of physical exercise in elderly people at highest risk of sarcopenia« . La qualité du reporting est bien pourrie, les critères d’inclusion flous, mais on peut comprendre qu’il s’agit d’une population de patients âgés fragiles, avec une sarcopénie mais en état clinique stable. Randomisation de 17 patients (groupe qui bénéfice de kinésithérapie seule) vs 17 patients (groupe qui bénéfice de kinésithérapie + FortiFit). Le FortiFit, est un complément alimentaire à base de protéines de lactosérum et de vitamines. Plusieurs tests standardisés sont passés à baseline et à trois mois. Furent enregistrés par impédancemétrie : la masse muculaire (kg), la masse maigre (kg), l’indice de masse maigre (kg/m²), la force musculaire au handgrip test, le Short Physical Performance Battery divisable lui-même en test d’équilibre, test de vitesse de marche et test de lever de chaise. Cela fait 8 critères de jugement potentiels, assez redondants. Trois méthodes statistiques principales sont aussi possibles pour les comparaisons du résultat à trois mois : test de Student sur les résultats à 3 mois, test de Student sur les changements (3 mois moins baseline), et modèle linéaire expliquant le résultat à 3 mois par le groupe de traitement et le résultat du test à baseline. Cela fait donc 24 analyses statistiques possibles. On ne sait pas trop quelles analyses étaient planifiées étant donné qu’on n’a pas accès au protocole.

Quel est le résultat ?

 Handgrip test Baseline Moyenne ± SEMHandgrip test 3 mois Moyenne ± SEMChangement moyen
Groupe contrôle23.73±2.06 kg22.18±2.19 kg-1.55
Groupe FortiFit22.51±2.35 kg24.54 ± 2.65 kg+2.03

Le test de Student des changements (3 mois moins baseline) est significatif (p=0.013), montrant que les sujets améliorent plus leur force musculaire avec le Handgrip test. Il y a une incohérence entre le texte et la figure 1 qui semble montrer seulement une augmentation moyenne de force de +1.6 kg dans le groupe FortiFit. Il y a peut-être un mélange accidentel de données entre les forces et masses musculaires (qui sont proches).

N’y a-t-il rien d’étonnant ? Pourquoi des patients, en état stable, perdent-ils 1.55 kg au handgrip test alors qu’ils bénéficient de kinésithérapie ? Pourquoi la différence entre les groupes à baseline est-elle dans le sens opposé de la différence finale ?

Pour commencer, considérons la distribution des changements moyens (ici, les changements moyens observés sont -1.55 kg et +2.03 kg), conditionnelle à une espérance nulle du changement dans chaque groupe. Sous hypothèse d’homoscédasticité, les deux changements moyens ont la même variance car les groupes sont de taille égale. Ils sont aussi indépendants. Par le théorème central limite, la distribution jointe de ces deux moyennes de changement est donc approximable à une distribution binormale avec une corrélation nulle, centrée autour du point (0,0).

La figure présente la distribution binormale ainsi que le seuil de significativité à 10% bilatéral (orange) et 5% bilatéral (rouge). Conditionnellement à un petit p < 0.05, on constate que la densité de proba bivariée est maximale pour des moyennes égales et de signe opposé. Présenté sous un autre angle, sachant que la différence est de 2 erreurs types entre les deux groupes, il est bien plus probable d’avoir -1 erreur type dans un groupe et +1 erreur type dans l’autre groupe, que +0 erreur type dans un groupe et +2 erreurs types dans l’autre et il est encore plus improbable d’avoir +1 erreur type dans un groupe et +3 erreurs types dans l’autre !

En conditionnant la distribution binormale à un petit p < 0.05, j’ai pu calculer numériquement une probabilité d’avoir un ratio des différences de moyenne compris entre -0.5 et -2. Cette probabilité est estimée à 56%. La probabilité d’un ratio négatif (c’est-à-dire, que les deux changements moyens sont de signe opposé) est estimée à 97.5%.

En réalisant la même considération en analyse tétravariée, on peut calculer que la densité de probabilité maximale est concentrée sur le scenario suivant

 Handgrip test Baseline Moyenne ± SEMHandgrip test 3 mois Moyenne ± SEMChangement moyen
Groupe contrôleµ + epsilonµ – epsilon-2×epsilon
Groupe FortiFitµ – epsilonµ + epsilon+2×epsilon

Où µ représente l’espérance commune aux quatre cases et epsilon est la valeur telle que la différence fournisse un p=0.05, c’est-à-dire la plus petite valeur qui conduise à un résultat statistiquement significatif.

Évidemment, le point de densité de probabilité le plus élevé reste infiniment improbable puisque les distributions sont continues. Les quatre écarts à la moyenne générale diffèreront donc plus ou moins, mais les grandes tendances devraient souvent apparaître.

Le cas observé semble donc bien typique, si on présente les différences par rapport à la moyenne générale (quatre cases moyennées) :

 Handgrip test Baseline Moyenne ± SEMHandgrip test 3 mois Moyenne ± SEMChangement moyen
Groupe contrôlem + 0.49m – 1.06 -1.55
Groupe FortiFitm – 0.73m + 1.30 +2.03

Un peu de pratique

Reprenons le cas d’école et appliquons la liste des critères :

  1. Échantillon de toute petite taille face à l’effet attendu, suggérant une puissance très faible. CHECK
  2. Estimation ponctuelle démesurée en comparaison à ce qui paraît plausible. DISCUTABLE
  3. Petit p à la limite de la significativité (p typiquement compris entre 0.01 et 0.05) CHECK
  4. Multiplicité des tests, apparente ou cachée CHECK. N=3 tests apparents (avec test complètement inapproprié) seulement, mais nombreux critères de jugements mesurés mais non présentés dans les résultats. Un des tests présentés n’était même pas dans ma liste des tests imaginables et cache un choix de seuillage sur une variable catégorielle ordinale (-> deux tests possibles).
  5. Autres tests répondant à la même question ne montrant pas plus de tendance à la significativité que ce qui est explicable par le hasard (environ un petit p sur 20 significatif, la moitié du temps dans le sens opposé à ce que veulent montrer les auteurs). FAILED. Sur les trois tests présentés, deux vont dans le même sens et le troisième n’est pas reproductible mais irait dans le même sens selon les auteurs. Ce FAIL pourrait être dû à la non présentation de nombreux tests cachés.
  6. Lorsqu’une différence d’évolution (p.e. Student inter-groupe sur différences intra-sujets) est analysée sur un paramètre dont la stabilité est attendue, apparition d’une dégradation importante dans le groupe contrôle et d’une amélioration de même amplitude dans le groupe expérimental et une différence à baseline qui va dans le sens opposé à la différence finale. CHECK pour le premier critère de jugement CHECK pour le second critère (masse musculaire) CHECK pour le troisième (sarcopénie binaire)

Sur la sémiologie de la multiplicité des tests :

  1. Écart au protocole, sur les analyses SUSPECT
  2. Critères de jugements présentés dans la partie méthodes mais pas dans les résultats CHECK
  3. Grande majorité (voire totalité) de petits p significatifs dans les résultats mais presque tous compris entre 0.01 et 0.05 CHECK
  4. Étrangeté des analyses qui « tirent dans les coins » FAILED
  5. Critère de jugement principal inattendu étant donné la population et l’intervention, voire disparition des critères de jugements attendus. FAILED

Conclusion

Avec un peu d’expérience, on peut identifier beaucoup de petits p bidons.

Hypothèse nulle et alternative

Vous avez peut-être entendu parler d’un résultat statistiquement significatif qui ne serait pas cliniquement significatif car correspondant à une différence trop faible, notamment à cause d’un échantillon « trop grand« . Le problème ne vient pas de la taille de l’échantillon mais du mauvais choix de l’hypothèse nulle. Si on veut prouver qu’un effet est cliniquement significatif, il faut que l’hypothèse alternative soit « cet effet dépasse le seuil de significativité clinique » et que l’hypothèse nulle en soit la négation, c’est-à-dire « cet effet est inférieur ou égal au seuil de significativité clinique ».

H0 : µ1-µ2 <= +clinthreshold

H1 : µ1-µ2 > +clinthreshold

Où clinthreshold représente le seuil de significativité clinique. On est habitué à ça pour les essais de non-infériorité, avec des seuils négatifs. Une analyse de supériorité devrait se faire de la même manière mais avec un seuil positif.

Évidemment, il ne faut pas se faire avoir par le syllogisme consistant à poser

H0 : µ1<=µ2

H1 : µ1 > µ2

puis rejeter H0 et conclure que la différence m1-m2 observée est la différence réelle et que, comme elle est supérieure au seuil de significativité clinique, on a prouvé qu’il existait un effet cliniquement significatif. En effet, dans le pire des cas, l’effet observé est égal au seuil de significativité clinique et l’affirmation selon laquelle l’effet réel est supérieur à ce seuil a un risque alpha unilatéral à 50%.

Pourquoi n’est-il pas pragmatique d’utiliser un seuil de significativité clinique ?

C’est trop subjectif. Comme tout ce qui est subjectif, les auteurs seraient tentés de tricher au maximum dessus (comme pour le delta du nombre de sujets nécessaires) et les reviewers pourraient toujours pinailler même si l’effet choisi par les auteurs est pertinent. Bref, ça pose des problèmes à tout le monde. Le zéro, par contre est objectif et consensuel même s’il est un des choix les moins pertinents qui soient.

Il y a d’autres problèmes, comme l’augmentation du nombre de sujets nécessaires (NSN), obligeant à tricher encore plus sur le calcul. La mascarade du NSN deviendrait encore plus évidente ; ce NSN est généralement calculé à l’envers, c’est-à-dire partant du nombre de sujets incluables pour en déduire les hypothèses nécessaires à générer ce nombre. Peut-être finirait-on par admettre qu’une étude ne répond pas de manière certaine à une question. Peut-être admettrait-on qu’il faut attendre la méta-analyse pour juger d’un effet et qu’un travail qui tente de répliquer un résultat a au moins autant de valeur sinon plus que l’article original, mais là, je rêve. Soyons réaliste, ça passerait mal.

Solution pragmatique : fournir un intervalle de confiance à 95% de l’effet et laisser le lecteur final de l’article déterminer lui-même si la borne basse dépasse ou non son seuil subjectif de significativité.

Poussons la réflexion

Vous l’aurez compris, ma formulation des hypothèses n’a qu’un intérêt conceptuel. Cela peut aider dans l’interprétation subjective de résultats et plus ou moins dans la construction méthodologique de projets de recherche, mais ne sera pas explicitement formulé.

Dans un essai clinique randomisé, on a généralement une hypothèse précise. On veut prouver la supériorité d’une intervention par rapport à une autre. Une vision « binaire » succès/échec est pertinente car la décision finale est binaire. En épidémiologie, on est parfois un peu plus neutre sur le sujet. Un cadre théorique à trois hypothèses paraît alors utilisable.

H0 : effet futile, c’est-à-dire -clinthresthold < effet < +clinthreshold

H1 : effet > clinthreshold

H2 : effet < clinthreshold

Si l’effet correspond à une différence entre deux groupes, alors on pourra reformuler :

H0 : équivalence

H1 : supériorité

H2 : infériorité

En analysant un intervalle de confiance, on pourra rejeter une ou plusieurs des trois hypothèses et fournir une conclusion plus ou moins fine.

PS : un intervalle de confiance avec un niveau de confiance différent est possible, bien sûr.

Six conclusions différentes sont possibles dont trois sont tranchées, deux sont un peu floues et une est très floue :

Conclusions tranchées:

Supériorité : on a rejeté H0 et H2 et on accepte H1

Équivalence : on a rejeté H1 et H2 et on accepte H0

Infériorité : on a rejeté H0 et H1 et on accepte H2

Conclusions un peu floues:

Non-infériorité : on a rejeté H2 mais H1 comme H0 restent compatibles avec les données

Non-supériorité : on a rejeté H1 mais H2 comme H0 restent compatibles avec les données

Conclusion très floue: on n’a rejeté aucune hypothèse. Toutes les hypothèses restent compatibles avec les données.

Dans ce cadre conceptuel, on pourrait conclure à un différence significativement futile si la conclusion est l’équivalence… quand bien même le zéro n’est pas contenu dans l’intervalle de confiance de la différence.

De mon point de vue, ce n’est pas forcément aux auteurs de l’article de déterminer les seuils de significativité clinique, mais plutôt au lecteur de l’article. Par ailleurs, on pourrait pousser le concept plus loin en distinguant les effets cliniquement significatifs mineurs des effets cliniquement significatifs majeurs, les deux n’impliquant pas forcément la même réaction. En bref, la meilleure manière de « lire » un intervalle de confiance, c’est de tenter de conclure séparément sur l’effet observé à chacune des deux bornes de l’intervalle de confiance et de se dire que la réalité est probablement quelque part entre les deux (en négligeant les problèmes de prior bayésien ; ce qui est acceptable quand l’étude est de suffisamment grande taille pour que l’information du prior soit négligeable face à l’information contenue dans les données de l’étude).

Délire unilatéraliste

Vous avez pu déjà voir la formulation unilatérale suivante

H0 : µ1 = µ2

H1 : µ1 > µ2

Cette formulation est incorrecte à moins qu’il y ait une preuve mathématique que la proposition µ1 < µ2 soit fausse de telle sorte que µ1=µ2 est mathématiquement équivalent à µ1<=µ2. C’est le cas pour le test de Fisher dans une ANOVA, mais c’est loin d’être une situation habituelle.

En effet, rejeter l’égalité ne prouve pas la supériorité. Si on construit une statistique égale à la valeur absolue de la différence entre m1 et m2, on peut arriver à la conclusion que µ1 n’est pas égal à µ2 et rejeter H0, mais ça ne prouvera nullement la supériorité !

Une formulation plus correcte serait :

H0 : µ1 <= µ2

H1 : µ1 > µ2

On peut s’apercevoir que la P-value n’est pas aisément calculable dans ce cas, parce que si µ1 est légèrement inférieur à µ2 ou très largement inférieur, alors la distribution de M1-M2 change beaucoup. Eh bien, dans ce contexte, on prendra la borne supérieure des P-values spécifiques d’une différence donnée µ1-µ2 négative ou nulle. Si la statistique est bien conçue, il est probable que le pire des cas (P-value la plus grande) est dans le scenario µ1=µ2. C’est peut-être pour ça que certains ont mal formulé les hypothèses. Néanmoins la bonne formulation des hypothèses évite la construction de statistiques buggées qui ne se comportent pas bien dans le cas où µ1<µ2.

Délire bilatéraliste

Vous avez certainement vu la formulation

H0 : effet = 0

H1 : effet ≠ 0

Cette fois-ci H1 est bien la négation de H0, mais deux problèmes apparaissent.

  1. L’hypothèse nulle est souvent invraisemblable. Par exemple, paraît-il possible que le traitement chirurgical et le traitement médical de la hernie discale aient exactement le même résultat fonctionnel moyen à 12 mois ? La différence peut-être minime, complètement négligeable, mais la nullité absolue de la différence paraît impossible.
  2. L’hypothèse alternative est pratiquement inutile car elle ne donne ni d’information sur l’amplitude, ni sur la direction de l’effet. Dire qu’une exposition a un effet n’a pas spécialement de pertinence si on ne précise pas si elle est bénéfique ou nocive…

Si les analyses statistiques sont simples, alors on peut s’aider de l’effet observé pour juger de la direction de l’effet. Ce n’est pas le cas avec certains modèles de survie. Un bon exemple est le test MaxCombo décrit par Theodore Karrison dans l’article intitulé « Versatile tests for comparing survival curves based on weighted log-rank statistics » publié dans « The Stata Journal » (2016, Vol 16, Number 3, pp. 678-690). Le test permet de rejeter l’hypothèse de superposition parfaite de deux courbes de survie (hypothèse nulle) mais n’aide pas à décider de laquelle est la meilleure. Si on conclut, par exemple, que le groupe pour lequel la plus grande médiane de survie est observée, a réellement une médiane de survie supérieure, on prend un risque alpha unilatéral pouvant atteindre 50% dans le pire des cas. De même pour la survie à 1 an, pour l’espérance de vie tronquée, pour l’espérance de vie globale et à peu près pour tout ce qu’on peut imaginer.

Conclusion

Hypothèse nulle ne devrait pas être synonyme d’absence totale d’effet mais devrait toujours être la négation de l’hypothèse alternative, cette dernière étant l’hypothèse que l’on souhaite prouver. Commencez toujours par formuler cette hypothèse alternative et vous produirez des hypothèses nulles pertinentes. On peut créer des cadres théoriques distinguant plus d’hypothèses, mais la dualité H0/H1 reste pertinente pour les essais cliniques.

Pourquoi un modèle multivarié ?

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

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

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

1. Analyse de causalité

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5. Recherche de « facteurs de risque »

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

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

À quoi cela sert-il ?

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

8. Imputation simple ou multiple

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

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

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

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

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

There will be legacy

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

Can a legacy programming language die?

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

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

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

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

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

Having many programming languages, create as many problems:

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

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

Very low level: assembly languages

Low level: C and C++

Intermediate with strong type system: C#, Java

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

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

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

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

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

Test de normalité

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

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

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

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

La normalité suppose aussi :

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

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

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

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

Hypothèse nulle H0 : le patient est immortel

Hypothèse nulle H1 : le patient est mortel

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Quelques néologismes

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

Petit péter, petit péteur

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

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

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

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

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

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

Le petit péteur a deux principales formes cliniques :

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

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

Chance alpha

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

Publicance

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

De manière plus formelle :

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

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

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

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

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

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

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

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

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

Facteur bayésien et puissance corrigée

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

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

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

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

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

Scenarii de simulation

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

Risque alpha bilatéral nominal : 5%

Outcome binaire : succès vs échec

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

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

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

Deux effets du traitement seront analysés :

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

Estimateurs comparés :

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

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

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

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

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

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

Fréquence de l’outcome

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

Gestion des défauts de convergence

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

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

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

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

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

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

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

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

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

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

Valeurs de départ et algorithme

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

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

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

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

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

Résultats modèle Gaussien

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

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

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

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

Résultats Wald sur identité-binomial

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

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

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

Résultats avec RV sur identité-binomial

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

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

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

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

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

Facteur bayésien

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

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

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

Défaut de convergence

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

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

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

Biais de l’estimation ponctuelle

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

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

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

Discussion

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

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

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

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

Combien de noeuds pour paralléliser mes simulations ?

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

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

Exemple : processeur 4 coeurs et 8 threads.

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

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

Limites de la règle n=threads

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

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

Pourquoi ?

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

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

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

Introduction

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

Résumé

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

Argumentaire détaillé

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

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

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

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

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

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

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

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

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

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

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