Auteur : Nicolas Stefaniak

Contributeur : Apolline Durtette

Niveau : Licence - master - doctorat - chercheurs confirmés

Prérequis

 install.packages("metafor")
 library("metafor")

Introduction

Nous avons vu jusqu’à présent ce qu’était une taille d’effet et comment les obtenir à partir des informations qui pouvaient être disponibles dans les articles scientifiques. L’étape suivante consiste à pouvoir combiner ces tailles d’effet. Cependant, nous sommes confrontés à deux problèmes :

  1. tout comme l’estimation des moyennes, des écart-types, de la valeur d’un statistique, quand on estime des tailles d’effet, la mesure que nous obtenons n’est pas une mesure parfaite. Il est donc nécessaire d’estimer l’erreur de mesure autour de ces tailles d’effet.

2)nous ne pouvons pas comparer directement les tailles d’effet en faisant simplement une moyenne car l’estimation de certaines tailles d’effet ont été obtenues sur de plus gros échantillons que d’autres, ce qui implique que les premières sont plus précises que les secondes.

Ansi, toute chose étant égale par ailleurs, une taille d’effet avec une plus petite estimation de son erreur sera plus précise. Les méta-analyses se doivent de prendre en compte cet aspect en accordant plus de poids à des tailles d’effets plus précises (Hedges & Olkin, 1985). un autre aspect qui va avoir un impact sur la précision est le design utilisé : avec une même taille d’échantillon, un échantillon apparié sera plus précis que deux échantillons indépendants, du moins si la corrélation entre les deux mesures n’est pas nulle (Borenstein et al., 2009). Á l’inverse, dans des modèles linéaires mixtes, le fait d’avoir des cluster augmente l’erreur-type, et donc diminue la précision, du moins si le coefficient intraclass n’est pas égal à 0.

Dans cette section, nous allons aborder la manière de prendre en compte les erreurs d’estimation dans les tailles d’effet lorsque nous faisons des méta-analyses en opposant deux logiques pour estimer ces erreurs d’estimation : les modèles à effets fixes et les modèles à effets aléatoires.

Les modèles fixes

Dans les modèles à effet fixe, nous faisons l’hypothèse que toutes les études inclues dans la méta-analyse partage une taille d’effet commune (Borenstein et al., 2009). Cela signifie que tous les facteurs qui ont un impact sur l’estimation de la taille d’effet sont identiques entre les études, ce qui implique que la taille d’effet est la même dans toutes les études, d’où l’appelation ‘fixe’. Cet effet fixe sera appelé \(\theta\).

Dans ce cas, les variations qui existent dans l’estimation de la taille d’effet s’explique simplement par une erreur aléatoire inhérente à chaque étude. Si chaque étude avait une taille d’échantillon énorme, l’erreur de mesure vaudrait 0. Il est évident que les tailles d’échantillons dans les études ne sont jamais suffisamment grandes pour pouvoir avoir une erreur d’estimation qui soit nulle.

Le graphique 1 est une représentation graphique de ce que nous venons de décrire. Dans ce graphique, les rectangle rouges représentent la taille d’effet réelle tandis que les point noirs repréentent la taille d’effet observée. Les tailles d’effet réelles sont parfaitement alignées sur l’axe vertical tandis que les points noirs se dispersent de manière aléatoire autour de ces points rouges. L’écart entre chaque point rouge et chaque point noir, représenté par les flèches, est l’erreur d’estimation de la taille d’effet dans notre modèle.

Figure 1. Représentation graphique de la taille d’effet réelle de 3 études (rectangle rouge) et de la taille d’effet observée pour ces études (points noirs)

Dans ce cas, l’idée est que, même si l’erreur se distribue de manière aléatoire, il est possible d’estimer la distribution d’échantillonnage de ces erreurs.

Le problème est que nous n’avons pas accès à la taille d’effet de la population, mais seulement à la taille d’effet observée dans chacune des étude. Il faut donc pouvoir estimer la taille d’effet de la population à partir de ces données. Pour obtenir l’estimation la plus précise possible de la taille d’effet, il faut calculer une moyenne pondérée de la taille d’effet ou le poids attribué à chaque étude doit être inversément proportionnel à l’erreur de mesure de la taille d’effet 1.

\[W_i = \frac{1}{V_{y_i}}\text{ (1)}\]

où \(V_yi\) est la variance la variance de chaque étude. On peut donc calculer la moyenne pondérée grâce à l’équation 2:

\[M = \frac{\displaystyle \sum_{i=1}^k W_i\times Y_i}{\displaystyle \sum_{i=1}^k W_i}\text{ (2)}\] où le \(W_i\times Y_i\) est le produit entre la taille d’effet et son poids, et \(W_i\) est le poids des études.

Cette taille d’effet globale est également associée à une erreur de mesure. La variance autour de la taille d’effet globale est obtenue par l’équation 3:

\[V_M = \frac{1}{\displaystyle \sum_{i=1}^k W_i}\text{ (3)}\]

l’erreur-type est donc obtenue en faisant la racine carrée de \(V_M\).

\[SE_M = \sqrt{V_M} \text{ (4)}\]

Á partir de là, il est très facile d’obtenir les limites supérieures et inférieures de l’intervalle de confiance autour de la taille d’effet réelle :

\[IC_{2.5} = M-1.96\times SE_M \text{ ET } IC_{97.5} = M+1.96\times SE_M\text{ (5)}\] On peut également tester la significativité de la taille d’effet, c’est-à-dire déterminer si la valeur de la taille d’effet est significativement différente de 0, grâce à un test Z :

\[Z =\frac{M}{SE_M}\text{ (6)}\]

La valeur de cette statistique se distribue en suivant une distribution normale.

Les modèles aléatoires

Dans le smodèles fixes, l’hypothèse qui est faites est que la taille d’effet est identique entre les études. Cependant, cette hypothèse est généralement peu plausible. Une hypothèse plus raisonnable est de considérer que les études ont suffisamment de choses en commun pour que cela fasse sens de les synthétiser, sans qu’on puisse pour autant faire l’hypopthèse que la taille d’effet est identique entre les études.

Les choix méthodologiques, les outils utilisés, les caractéristiques de la population testée … peuvent varier d’une étude à l’autre, et ces différences peuvent consécutivement impacter l’estimation de la taille d’effet. Par ailleurs, il peut y avoir toute une série de caractéstiques environnementales qui peuvent avoir un impact sur la taille d’effet.

Par exemple, si on s’intéresse à la pertinence d’une nouvelle approche pédagogique, d’une classe à l’autre, d’une école à l’autre, il peut y avoir des ressources qui varient et il est tout à fait vraisemblable que ces ressources n’aient pas été examinées d’une étude à l’autre.

Pour prendre en compte l’ensemble des facteurs qui peuvent impacter l’estimation de la taille d’effet, on peut réaliser une méta-analyse en choisissant un effet aléatoire. L’hypothèse sous-jacente est que la taille d’effet est distribuée en suivante un distribution normale entre les études.

Le graphique 2 est une représentation graphique de ce que nous venons de décrire. Dans ce graphique, les rectangle rouges représentent la moyenne des taille d’effet réelle pour l’ensemble des études, les triangles verts représentent la taille d’effet réelle pour chacune des études et les points noirs représentent les tailles d’effet observées.

Figure 2. Représentation graphique de la taille d’effet réelle de 3 études (triangles verts) de la taille d’effet observée pour ces études (points noirs) et de la taille d’effet moyenne de chacune des tailles d’effet individuelles

Cette Figure illustre le fait que la distance entre la taille d’effet commune et les tailles d’effet observées se répartit en deux parties : les variations de taille d’effet entre les études \(\zeta_i\) représenté par les flèches noires, et les erreurs d’échantillonnage, représentées par les flèches bleues.

D’un point de vue mathématique, on peut donc représenter cette conception par l’équation

\[Y_i =\mu + \zeta_i+\epsilon_i\text{ (7)}\]

Nous devons donc non seulement évaluer la variance d’erreur, comme dans le modèle fixe, mais également la variance dans la distribution des tailles d’effet réelle.

La distance entre la taille d’effet moyenne commune \(\mu\) et la taille d’effet réelle spécifique à chaque étude dépend de l’écart-type de la distribution des tailles d’effet est appelée \(\tau\), et sa variance est \(\tau^2\). Cette valeur s’applique pour l’ensemble des études inclues dans la méta-analyse.

La distance entre chacune des tailles d’effet réelle \(\theta_i\) et la taille d’effet observée \(Y_i\) dépend de la variance \(V_i\) qui va varier être spécifique à chacune des études.

Á partir des tailles d’effet observée, nous voulons estimer la taille d’effet commune \(\mu\). Tout comme dans les modèles à effet fixe, nous commençons par calculer une moyenne pondérée des tailles d’effet, où la pondération est l’inverse de la variance de chaque étude.

Pour calculer les variances dans les modèles à effet aléatoire, on va appliquer le principe de décomposition de la variance : la variance totale est la somme de la variance spécifique à chaque étude et de la variance des effets réels \(\tau^2\).

Pour estimer \(\tau^2\), nous devons utiliser la méthode des moments (ou la méthode DerSimonian) (Borenstein et al., 2009) présentée dans l’équation 8.

\[\tau^2 =\frac{Q-df}{c} \text{ (8)}\]

où Q vaut :

\[ Q =\displaystyle \sum_{i=1}^k W_i\times Y_i^2 - \frac{(\displaystyle \sum_{i=1}^k W_i\times Y_i)^2}{\displaystyle \sum_{i=1}^k W_i} \text{ (9)}\]

Les degrés de liberté (df) correspondent au nombre d’études moins 1 et C vaut :

\[ C =\sum W_i\ - \frac{\sum W_i^2}{ \sum W_i} \text{ (10)}\]

Pour estimer la taille d’effet commune, on va appliquer, comme dans les modèles fixes, une pondération à la moyenne 11. L’astérisque signifie que nous sommes dans un modèle aléatoire.

\[W^*_i=\frac{1}{V^*_{Yi}} \text{ (11)}\]

où \(V^*_{Yi}\) est la somme de la variance intra-étude et entre les études.

\[V^*_{Yi}= V_{Yi}+\tau^2 \text{ (12)}\]

On obtient donc la moyenne pondérée par

\[M^* = \frac{\displaystyle \sum_{i=1}^k W^*_i\times Y_i}{\displaystyle \sum_{i=1}^k W^*_i}\text{ (13)}\]

La variance sera estimée par la récuproque de la somme des pondérations.

\[V_{M^*} = \frac{1}{\displaystyle \sum_{i=1}^k W^*_i}\text{ (14)}\]

l’erreur-type est donc obtenue en faisant la racine carrée de \(V_{M^*}\).

\[SE_{M^*} = \sqrt{V_{M^*}} \text{ (15)}\]

Á partir de là, il est très facile d’obtenir les limites supérieures et inférieures de l’intervalle de confiance autour de la taille d’effet réelle :

\[IC_{2.5} = M^*-1.96\times SE_{M^*} \text{ ET } IC_{97.5} = M^*+1.96\times SE_{M^*}\text{ (16)}\] On peut également tester la significativité de la taille d’effet, c’est-à-dire déterminer si la valeur de la taille d’effet est significativement différente de 0, grâce à un test Z :

\[Z^* =\frac{M^*}{SE_{M^*}}\text{ (17)}\]

La valeur de cette statistique se distribue en suivant une distribution normale.

Faut-il choisir un effet fixe ou un effet aléatoire ?

Une différence concrète quand on choisit entre un modèle fixe et un modèle aléatoire est que, lorsque nous faisons un modèle à effet fixe, nous pouvons largement ignorer l’effet des plus petites études étant donné que nous allons avoir une meilleure appréciation de l’information dans les études avec de plus gros échantillons. Dans les modèles à effet aléatoire, nous ne souhaitons pas représenter un effet mais la moyenne de la distribution des vrais effets. Il faut donc que toutes les tailles d’effet soient représentées. Cela signifie qu’on ne peut ignorer une petite étude en lui attribuant un tout petit poids, puisque même si cette taille d’effet est imprécise, c’est l’estimation d’une taille d’effet qu’aucune autre étude n’a estimé. Á l’inverse, on ne peut pas donner un poids prépondérant à une étude avec un très gros échantillon. Cela signifie donc que dans les modèle aléatoires, les poids entre les études seront beaucoup plus restreints que dans un modèle à effet fixes. Lorsque nous réaliserons les analyses avec le logiciel, je vous encourage à comparer le poids des études selon que le modèle est à effet fixe ou à effet aléatoire.

En d’autres termes, il y aura un équilibre plus marqué entre les études dans un modèle à effet aléatoire. Donc, les grandes études vont perdre de leur influence alors que les petites vont en gagner. Cela est donc d’autant plus vraies pour les études qui ont des valeurs extrêmes.

Une autre conséquence est que la variance d’erreur sera systématique plus large pour les modèles à effets aléatoires que pour les modèles à effets fixes. Ceci s’explique assez aisément : imaginez que la taille de votre échantillon soit extrêmement grande, la variance d’erreur va tendre vers 0. Dans un modèle à effet fixe, cette variance de la taille d’effet sera donc de 0. En revanche, dans les modèles à effet aléatoire, si la variance pour chaque étude va tendre vers 0, la variation entre les études continuera à persister. Pour que la variance entre les études soit proche de 0, il faudrait que le nombre d’études inclues dans la méta-analyse tende vers l’infini.

Aspects mathématiques expliquant en quoi l’erreur de mesure est plus large pour un modèle aléatoire que pour un modèle fixe

Lorsque nous réalisons un modèle à effet fixe, l’erreur-type vaut :

\[SE_M=\sqrt{\frac{\sigma^2}{k \times n}}\] où k est le nombre d’études et n est le nombres de participants. Ainsi, quand n grandit, \(SE_M\) diminue pour tendre vers 0 quand n tend vers l’infini.

Dans un modèle à effet aléatoire, l’erreur-type est obtenue de la manière suivante :

\[SE_M=\sqrt{\frac{\sigma^2}{k \times n}+\frac{\tau^2}{k}} \]

Ainsi, il faut non seulement que n tend vers l’infini mais aussi que le nombre d’études tend vers l’infini puisque, si la première partie est identique à ce qui se passe dans les modèles à effet fixe, la seconde partie réflète la variance entre les études.

Les deux types de modèles se distinguent également sur l’hypothèse nulle où l’hypothèse nulle des modèles à effets fixes est que la taille d’effet moyenne est égale à 0, alors que dans les modèles à effets aléatoires, l’hypothèse nulle est que la moyenne des tailles d’effet vaut 0, ce qui signifie que les vraies tailles d’effet des études individuelles peuvent être différentes de 0.

Alors comment choisir entre effets fixes ou effets aléatoires ? On peut choisir des modèles à effets fixes si nous avons la conviction que toutes les études inclues dans la méta-analyse sont identiques et que notre objectif est de calculer la taille d’effet commune que nous ne souhaitons pas généraliser à une autre population que celle investiguée. Si les études ont été réalisées par des chercheurs qui ont mené des études de manière indépendante, il est peu vraisemblable que toutes les études soient fonctionnellement équivalentes (Borenstein et al., 2009). En d’autres termes, si on considère que les tailles d’effet comme ayant été échantillonnées à partir d’une distribution de taille d’effet, les modèles aléatoires doivent être préférés.

De plus, si la variance entre les études et importantes, et significative, les modèles à effets fixes sont inappropriés.

En psychologie, il semble y avoir peu de cas où il serait raisonnable d’utiliser des modèles à effets fixes plutôt que des modèles à effet aléatoire. En effet, même si nous avons été attentif à choisir des études parfaitement homogènes en utilisant la méthodologie PICO, c’est peu vraisemblables que les études soient parfaitement homogènes. Il est donc concentionnel de toujours utiliser des modèles aléatoires (Harrer et al., 2021). Les effets fixes n’auraient de sens que dans les situations où il s’agit d’études reproduites à l’identique, ce qui est rarement le caS.

Si les modèles à effets aléatoires semblent naturels en psychologie, le fait que les petites études reçoivent plus de poids dans des modèles à effets aléatoires que dans des modèles à effets fixes peut être problématique car elles sont plus susceptibles d’être empreintes de biais (Borenstein et al., 2009; Schwarzer et al., 2015). C’est pourquoi certains avancent qu’il est parfois préférable d’utilser des modèles à effets fixes (Furukawa et al., 2003; Poole & Greenland, 1999) ou des modèles basés sur la méthode des moindres carré pondéré sans restriction (en anglais : unrestricted weighted least squares - UWLS)(Stanley et al., 2022).

Réalisation d’une méta-analyse à l’aide de R.

Les données pour les exemples de travail sont présentées dans le fichier excel ci-joint.

{{% attachments style="grey" pattern=".*\.xlsx$" %}}

Dans R, il existe plusieurs packages en mesure de réaliser des méta-analyses. Parmi ces packages, il y a :

Ces trois packages sont très bien faits et fournissent des résultats sensiblement identiques. Le choix pour le package ‘metafor’ a été fait pour la présentation principale des analyses dans le cadre de ce document. La raison de ce choix est que ‘metafor’ est le plus complet des packages, et en particulier que ‘meta’. Par ailleurs, ‘dmetar’ n’est pas un package officiel de R. Il n’a donc pas de documentation pour chaque ligne de commande comme c’est le cas pour les packages officiels de R.

Pour les personnes qui estimeraient que l’utilisation de ces fonctions sont compliquées, vous pouvez envisager les alternatives qui sont présentées.

Pour la mise en application, il est nécessaire que les packages ‘readxl’ et ‘metafor’ soient installés et chargés.

install.packages("readxl") # à ne faire que si ce package n'est pas déjà installé
install.packages("metafor")# à ne faire que si ce package n'est pas déjà installé

require(readxl)
require(metafor)

Ensuite, on choisit le répertoire de travail dans lequel nous souhaitons travailler. Ce répertoire doit nécessaire être celui dans lequel se situe le fichier de données.

getwd() # obtenir le répertoire de travail
setwd("C:/Mon_Répertoire")

Astuce

Pour éviter de devoir spécifier le répertoire de travail, au lieu de commencer en créant un script, vous pouvez commencer en créant un nouveau projet. Vous enregistrez le projet dans le dossier dans lequel se trouve les données et Rstudio fixera par défaut le répertoire de travail dans le répertoire sélectionné.

Les méta-analyses sur des différences de moyennes standardisées

On importe les données avec la fonction read_excel. Les données pour cet exercice sont présentées dans la feuille “ex1”. L’exercice est tiré de Borenstein et al. (Borenstein et al., 2009). Cet exemple permet d’illustrer une méta-analyse sur des différences de moyennes. Borestein et al. ne fournissent pas plus d’informations sur ces données.

dat <- read_excel("Exercices.xlsx", sheet="ex1")

Les données sont présentées dans le jeu de données ci-dessous.

Exemple issu de Bornstein et al. 2009

study

mean1

sd1

n1

mean2

sd2

n2

Carroll

94

22

60

92

20

60

Grant

98

21

65

92

22

65

Peck

98

28

40

88

26

40

Donat

94

19

200

82

17

200

Stewart

98

21

50

88

22

45

Young

96

21

85

92

22

85

Dans ce jeu de données, nous avons dans la première colonne le nom des études, ensuite les moyennes, écart-type, et taille d’échantillon pour chacun des groupes. Les tailles d’effet ne sont pas disponibles. La première étape va consister à ajouter les tailles d’effet à notre jeu de données. Pour obtenir les tailles d’effet, nous pouvons utiliser la fonction esclc. Pour obtenir les tailles d’effet, un certain nombre d’arguments vont devoir être précisés.

  • measure : il s’agit de la mesure que nous désirons. Pour des différences de moyennes, nous pouvons vouloir comparer :
  1. les différences de moyennes brutes - “MD”,
  2. les différences de moyennes standardisées (G de Hedge) - “SMD”,
  3. les différences de moyennes standardisées pour des variances hétérogènes - “SMDH”,
  4. les différences de moyennes standardisées avec le second groupe comme groupe de référence \(\Delta\) de Glass - SMD1
  5. les différences de moyennes standardisées avec le second groupe comme groupe de référence pour des variances hétérogènes - SMD1
  • m1i : la variable dans le jeu de données où se situe la moyenne du groupe 1
  • sd1i : la variable dans le jeu de données où se situel’écart-type du groupe 1
  • n1i : la variable dans le jeu de données où se situe la taille de l’échantillon du groupe 1
  • m2i : la variable dans le jeu de données où se situe la moyenne du groupe 2
  • sd2i : la variable dans le jeu de données où se situe l’étcart-type du groupe 2
  • n2i : la variable dans le jeu de données où se situe la taille de l’échantillon du groupe 2
  • slab : la variable dans le jeu de données où se situe l’identifiant des études (cet argument est accessoire)
  • data : le nom du jeu de données

Dans notre exemple, nous souhaitons réaliser une méta-analyse sur des différences de moyennes standardisées. Remarquez qu’après l’exécution de la fonction, il y a 2 nouvelles variables, yi et vi, qui sont respectivement les tailles d’effets et les erreur-types pour ces tailles d’effet.

dat <- escalc(measure= "SMD", m1i=mean1, sd1i=sd1, n1i=n1, m2i=mean2, sd2i=sd2, n2i=n2,
              slab=study, data=dat)

Remarquez à présent que les tailles d’effet et les erreurs-types ont été intégrées au jeu de données

Exemple issu de Bornstein et al. 2009 - La variable yi est la taille d'effet et la variable vi est l'erreur-type de la taille d'effet

study

mean1

sd1

n1

mean2

sd2

n2

yi

vi

Carroll

94

22

60

92

20

60

0.09452416

0.03337056

Grant

98

21

65

92

22

65

0.27735587

0.03106510

Peck

98

28

40

88

26

40

0.36654443

0.05083972

Donat

94

19

200

82

17

200

0.66438497

0.01055176

Stewart

98

21

50

88

22

45

0.46180628

0.04334467

Young

96

21

85

92

22

85

0.18516444

0.02363025

Nous voulons à présent réaliser la méta-analyse en tant que telle. On réalise cette méta-analyse avec la fonction rma. Cette fonction requiert 4 arguments :

  1. yi : la colonne avec la taille d’effet
  2. vi : la colonne avec la variance de la taille d’effet
  3. data : le jeu de données où on trouve les variables précisée pour yi et vi
  4. method : la méthode à devoir utiliser. Plusieurs possibilités existent :
    1. method = “EE” : les modèles à effets égaux ;
    2. method = “FE” : les modèles à effets fixes ;
    3. method = “DL” : modèle à effets aléatoires en utilisant la méthode DerSimonian-Laird pour réaliser l’estimation ;
    4. method = “HE” : modèle à effets aléatoires en utilisant la méthode de Hedges pour réaliser l’estimation ;
    5. method = “HSK” : modèle à effets aléatoires en utilisant lla méthode de Hunter-Schmidt spécifiquement adaptée pour les petits échantillons pour réaliser l’estimation ;
    6. method = “SJ” : modèle à effets aléatoires en utilisant la méthode de Sidik-Jonkman pour réaliser l’estimation ;
    7. method = “ML” : modèle à effets aléatoires en utilisant le maximum de vraisemblance pour réaliser l’estimation ;
    8. method = “REML” : modèle à effets aléatoires en utilisant le maximum de vraisemblance restreint pour réaliser l’estimation ;
    9. method = “EB” : modèle à effets aléatoires en utilisant l’estimateur empirique de Bayes ;
    10. method = “PM” : modèle à effets aléatoires en utilisant la méthode Paule-Mandel pour réaliser l’estimation ;
    11. method = “GENQ” : modèle à effets aléatoires en utilisant la statistique Q généralisée pour réaliser l’estimation ;
    12. method = “PMM” : modèle à effets aléatoires en utilisant la méthode de Paule-Mandel basée sur la médiane non biaisée pour réaliser l’estimation ;
    13. method = “GENQM” : modèle à effets aléatoires en utilisant la médiane non biaisée de la statistique Q généralisée pour réaliser l’estimation ;
    14. method = “HS” : modèle à effets aléatoires en utilisant de Hunter-Schmidt pour réaliser l’estimation.

Le méthdodes c à m sont des algorithmes qui réaliser des modèles à effets aléatoires. Parmi ces méthodes, la méthode de Paule-Mandel est la meilleur pour les variables dichotomiques et continues alors que le maximum de vraisemblance restreint est la meilleure pour les variable quantitatives. La méthode la plus adaptée pour obtenir l’intervalle de confiance est la méthode “GENQ” (Veroniki et al., 2015).

Il faut noter que les modèles à effets égaux et les modèles à effets fixes s’appuient sur les mêmes algorithmes. Ce qui les distingue, c’est que l’hypothèse qui sous-tend les deux modèles. Pour les modèles à effets égaux fait l’hypothèse que l’effet réel est homogène. Si nous rejetons l’hypothèse d’homogénéité, alors on doit rejeter le modèle dans son ensemble. Dans les modèles à effets fixes on ne fait pas l’hypothèse d’homogénéité mais on in interprète l’estimation comme une moyenne pondérée de l’effet réel des études Laird & Mosteller (1990).

Nous allons commencer par réaliser un modèle à effet fixe.

model_EE <- rma(yi = yi, # la colonne dans le jeu de données où on trouve la taille d'effet 
                vi = vi, # la colonne où on trouve l'erreur de mesure  
                data=dat, # le nom du jeu de données 
                method="EE", # la méthode permet de déterminer si c'est un effet fixe ou aléatoire. L'option "EE" renvoie à un effet fixe 
                digits=2) # le nombre de décimale à imprimer dans la sortie de résultats
summary(model_EE)
## 
## Equal-Effects Model (k = 6)
## 
##   logLik  deviance       AIC       BIC      AICc   
##    -0.83     11.91      3.66      3.45      4.66   
## 
## I^2 (total heterogeneity / total variability):   58.03%
## H^2 (total variability / sampling variability):  2.38
## 
## Test for Heterogeneity:
## Q(df = 5) = 11.91, p-val = 0.04
## 
## Model Results:
## 
## estimate    se  zval  pval  ci.lb  ci.ub      
##     0.42  0.06  6.46  <.01   0.29   0.54  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dans cette sortie de résultats, nous pouvons déjà focaliser notre attention sur l’estimate qui est la taille d’effet moyenne. Elle vaut 0.42 et la statistique Z qui lui est associée vaut 6.46, ce qui est significatif. Notre vraie taille d’effet se situe quelque part entre 0.29 et 0.54.

Remarquez que la pondération dont nous avons parlé dans la partie théorique est disponible de la manière suivante :

weights(model_EE)
##   Carroll     Grant      Peck     Donat   Stewart     Young 
## 12.383524 13.302553  8.128392 39.163627  9.533933 17.487971

Habituellement, on représente les méta-analyses sous la forme d’un graphique en forêt (forest plot en anglais). On obtient ce type de graphique avec la fonction forest Pour ce forest plot, nous souhaiterions afficher le poids des études dans le graphique.Nous venons de voir que les pondérations sont des nombres sans le symbole “%”, et qu’il ne sont pas arrondi. On va commencer par préparer la présentation des poids pour avoir des nombres arrondis suivi du symbole “%”. Pour cela, on va utiliser la fonction round pour faire l’arrondi et la fonction paste0 pour concatener les informations.

dat$weights <- paste0(round(weights(model_EE)), "%") # weights in % (rounded)

dat$weights 
## [1] "12%" "13%" "8%"  "39%" "10%" "17%"

Pour réaliser un forest plot, plusieurs étape sont nécessaires. Chacune des étapes va être commentée dans le code.

 # paramètres des graphiques. Les 4 valeurs indiquent les marges à appliquer au grpahique 
 # Dans l'ordre le bas, gauche, haut et droite

par(mar=c(4,4,2,2)) 

forest(model_EE, # le graphique en forêt doit être fait sur la sortie de la méta-analyse 
       xlim=c(-2,3), # les limites horizontales de la région du graphique. 0,0 est le centre du graphique
       header=TRUE, # faut-il ou non afficher le nom des colonnes au graphique
       top=2, # le nombre de ligne à laisser blanche en haut du graphique (pour pouvoir y ajouter des choses ensuite)
       mlab="Summary", # le nom à afficher pour la la ligne de synthèse des tailles d'effet 
       ilab=dat$weights, # les colonnes qu'on veut rajouter dans le graphique et qui ne sont pas présentes par défaut. On rajoute le poids des études ici. 
       ilab.xpos=1.6, # la position où il faut rajouter les colonnes supplémentaires. Pour rappel 0 est le centre du graphique
       ilab.pos=2) # l'alignement souhaité. Si on ne le précise pas, c'est centré. 2 = aligné à gauche 
text(x=1.6, y= -1, "100%", pos=2) # on rajoute en position x = 1.6, y = -1, la valeur 100% 
text(x= 1.6,  y= 8, "Weight", pos=2, font=2) # on rajoute le mot "weight" en position x=1.6 et y = 8

Astuce

Pour exporter ce graphique avec une qualité d’image suffisante pour un article scientifique, on peut utiliser la fonction png. Le premier argument est le nom du fichier dans lequel l’image sera enregistrée. Dans cette fonction, on peut préciser la hauteur et la largeur de l’image, ainsi que la résolution. Pour un article scientifique, il faut une résolution minimale de 300 dpi (dot per inch - point par pouce). En l’occurrence, 400 dpi a été choisi.

png(file='forestplot.png',width= 3000, height=1500, res = 400 ) # Open PNG device with specific file name
par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(model_EE)), "%") # weights in % (rounded)

forest(model_EE, xlim=c(-2,3), header=TRUE, top=2, mlab="Summary",
       ilab=dat$weights, ilab.xpos=1.6, ilab.pos=2)
text(1.6, -1, "100%", pos=2)
text(1.6,  8, "Weight", pos=2, font=2)
dev.off()

Exercice 1

Le Syndrome d’Apnées Obstructives du Sommeil (SAOS) occasionne, par ses mécanismes de pression thoracique et de désaturation en oxygène, de plus grands efforts de la part du cœur pour maintenir les fonctions vitales lors des différentes phases du sommeil. Cette charge mécanique entraine une plus grande fragilité de cet organe. À long terme, on peut s’attendre à une augmentation significative des risques cardiovasculaires chez les personnes souffrant de SAOS.

Les auteurs se sont demandés si la propagation des marqueurs de l’inflammation immunitaire au sein de l’organisme était un facteur de risque d’un trouble cardiovasculaire. Les auteurs se demandent si la réaction inflammatoire était significativement réduite chez les personnes qui ont bénéficié d’un traitement par Pression Positive Continue (PPC) puisque cette thérapeutique réduit le risque de trouble cardiovasculaire.

Les données sont disponibles sur la feuille ‘ex4’. Importez les données, calculez les tailles d’effet, réaliser la méta-analyse et présentez les résultats sous la forme d’un forest plot (en ayant exporté le fichier png). Faites cet exercice en réalisant l’analyse avec un effet fixe, ensuite avec un effet aléatoire et comparer les résultats.

Cliquez pour la solution 
dat <- read_excel("Exercices.xlsx", sheet="ex4")

Les données sont présentées ci-dessous.

Données SAOS - différences de moyennes standardisées

study

inf_index_baseline

sd_baseline

n_baseline

inf_index_ttmt

sd_ttmt

n_ttmt

Blue et al., 2020

80.2

5.21

136

83.0

3.85

136

Gulia et al., 2018

79.5

7.32

103

74.6

6.24

100

Jackson et al., 2019

87.6

2.01

128

85.7

2.05

126

Jonas et al., 2017

85.1

5.11

34

81.0

4.16

34

Leng et al., 2022

82.0

6.58

92

79.2

6.50

87

Schroeder et al., 2019

70.6

2.30

76

72.4

6.49

74

Thompson et al., 2022

73.9

7.69

120

70.3

5.60

119

Wolf et al., 2019

70.1

2.35

100

72.5

2.14

98

Zhunik et al., 2017

76.9

3.08

87

73.2

4.02

87

##### > Étape n° 4 : Calcul des variances et des tailles d'effet 
dat <- escalc("SMD", m1i=inf_index_baseline, sd1i=sd_baseline, n1i=n_baseline, 
              m2i=inf_index_ttmt, sd2i=sd_ttmt, n2i=n_ttmt,
              slab=study, data=dat)
dat
## 
##                    study inf_index_baseline sd_baseline n_baseline 
## 1      Blue et al., 2020               80.2        5.21        136 
## 2     Gulia et al., 2018               79.5        7.32        103 
## 3   Jackson et al., 2019               87.6        2.01        128 
## 4     Jonas et al., 2017               85.1        5.11         34 
## 5      Leng et al., 2022               82.0        6.58         92 
## 6 Schroeder et al., 2019               70.6        2.30         76 
## 7  Thompson et al., 2022               73.9        7.69        120 
## 8      Wolf et al., 2019               70.1        2.35        100 
## 9    Zhunik et al., 2017               76.9        3.08         87 
##   inf_index_ttmt sd_ttmt n_ttmt      yi     vi 
## 1           83.0    3.85    136 -0.6096 0.0154 
## 2           74.6    6.24    100  0.7169 0.0210 
## 3           85.7    2.05    126  0.9332 0.0175 
## 4           81.0    4.16     34  0.8699 0.0644 
## 5           79.2    6.50     87  0.4262 0.0229 
## 6           72.4    6.49     74 -0.3698 0.0271 
## 7           70.3    5.60    119  0.5331 0.0173 
## 8           72.5    2.14     98 -1.0633 0.0231 
## 9           73.2    4.02     87  1.0287 0.0260
##### > Étape n° 5.a : réalisation du modèle à effet fixe (Equal Effect model) 
model_EE <- rma(yi, vi, data=dat, method="EE", digits=2)
model_EE
## 
## Equal-Effects Model (k = 9)
## 
## I^2 (total heterogeneity / total variability):   96.17%
## H^2 (total variability / sampling variability):  26.13
## 
## Test for Heterogeneity:
## Q(df = 8) = 209.05, p-val < .01
## 
## Model Results:
## 
## estimate    se  zval  pval  ci.lb  ci.ub      
##     0.22  0.05  4.50  <.01   0.13   0.32  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Traduction des résultats en un "forestplot", incluant le poids ("weights")
## de chaque article dans les résultats
par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(model_EE)), "%") # weights en % (arrondis)

forest(model_EE, xlim=c(-4,4), header=TRUE, top=2, mlab="Summary",
       ilab=dat$weights, ilab.xpos=2.2, ilab.pos=2)
text(2.2, -1, "100%", pos=2)
text(2.2,  11, "Weight", pos=2, font=2)

Nous pouvons à présent réaliser un modèle avec un effet aléatoire.

##### > Étape n° 4 : Calcul des variances et des tailles d'effet 
dat <- escalc("SMD", m1i=inf_index_baseline, sd1i=sd_baseline, n1i=n_baseline, 
              m2i=inf_index_ttmt, sd2i=sd_ttmt, n2i=n_ttmt,
              slab=study, data=dat)
dat
## 
##                    study inf_index_baseline sd_baseline n_baseline 
## 1      Blue et al., 2020               80.2        5.21        136 
## 2     Gulia et al., 2018               79.5        7.32        103 
## 3   Jackson et al., 2019               87.6        2.01        128 
## 4     Jonas et al., 2017               85.1        5.11         34 
## 5      Leng et al., 2022               82.0        6.58         92 
## 6 Schroeder et al., 2019               70.6        2.30         76 
## 7  Thompson et al., 2022               73.9        7.69        120 
## 8      Wolf et al., 2019               70.1        2.35        100 
## 9    Zhunik et al., 2017               76.9        3.08         87 
##   inf_index_ttmt sd_ttmt n_ttmt      yi     vi weights 
## 1           83.0    3.85    136 -0.6096 0.0154     16% 
## 2           74.6    6.24    100  0.7169 0.0210     12% 
## 3           85.7    2.05    126  0.9332 0.0175     14% 
## 4           81.0    4.16     34  0.8699 0.0644      4% 
## 5           79.2    6.50     87  0.4262 0.0229     11% 
## 6           72.4    6.49     74 -0.3698 0.0271      9% 
## 7           70.3    5.60    119  0.5331 0.0173     14% 
## 8           72.5    2.14     98 -1.0633 0.0231     11% 
## 9           73.2    4.02     87  1.0287 0.0260      9%
##### > Étape n° 5.a : réalisation du modèle à effet fixe (Equal Effect model) 
model_EE <- rma(yi, vi, data=dat, method="PM", digits=2)
model_EE
## 
## Random-Effects Model (k = 9; tau^2 estimator: PM)
## 
## tau^2 (estimated amount of total heterogeneity): 0.55 (SE = 0.29)
## tau (square root of estimated tau^2 value):      0.74
## I^2 (total heterogeneity / total variability):   96.10%
## H^2 (total variability / sampling variability):  25.67
## 
## Test for Heterogeneity:
## Q(df = 8) = 209.05, p-val < .01
## 
## Model Results:
## 
## estimate    se  zval  pval  ci.lb  ci.ub    
##     0.27  0.25  1.06  0.29  -0.23   0.77    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Traduction des résultats en un "forestplot", incluant le poids ("weights")
## de chaque article dans les résultats
par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(model_EE)), "%") # weights en % (arrondis)

forest(model_EE, xlim=c(-4,4), header=TRUE, top=2, mlab="Summary",
       ilab=dat$weights, ilab.xpos=2.2, ilab.pos=2)
text(2.2, -1, "100%", pos=2)
text(2.2,  11, "Weight", pos=2, font=2)

On observe que le poids de l’étude de Blue était de 16% dans un modèle à effets fixes et il est de 11% dans le modèle à effets aléatoires. Á l’inverse, l’étude de Jonas avait un poids de 4% dans le modèle à effets fixes alors que son poids est de 10% dans ke modèle à effets aléatoires. Ces résultats permettent d’identifier la manière dont les modèles à effets aléatoires rééquilibrent le poids accordé à chacune des études.

Les méta-analyses sur des rapports de cotes

Des chercheurs se sont intéressés à l’efficacité du vaccin BCG contre la tuberculose. Trois études ont été inclues. Les données peuvent être chargées dans la mémoire de R en utilisant la commande

data(dat.bcg)

Les données sont présentées ci-dessous.

Données de l'efficacité du vaccin BCG contre la tuberculose

trial

author

year

tpos

tneg

cpos

cneg

ablat

alloc

1

Aronson

1,948

4

119

11

128

44

random

2

Ferguson & Simes

1,949

6

300

29

274

55

random

3

Rosenthal et al

1,960

3

228

11

209

42

random

4

Hart & Sutherland

1,977

62

13,536

248

12,619

52

random

5

Frimodt-Moller et al

1,973

33

5,036

47

5,761

13

alternate

6

Stein & Aronson

1,953

180

1,361

372

1,079

44

alternate

7

Vandiviere et al

1,973

8

2,537

10

619

19

random

8

TPT Madras

1,980

505

87,886

499

87,892

13

random

9

Coetzee & Berjak

1,968

29

7,470

45

7,232

27

random

10

Rosenthal et al

1,961

17

1,699

65

1,600

42

systematic

11

Comstock et al

1,974

186

50,448

141

27,197

18

systematic

12

Comstock & Webster

1,969

5

2,493

3

2,338

33

systematic

13

Comstock et al

1,976

27

16,886

29

17,825

33

systematic

Dans ce tableau, les variables correspondent à :

  • trial : le numéro de l’essai
  • author : le auteurs de l’étude
  • year : l’année de réalisation de l’étude
  • tpos : le nombre de personnes positives à tuberculose chez les personnes vaccinées
  • tneg : le nombre de personnes négatives à tuberculose chez les personnes vaccinées
  • cpos : le nombre de personnes positives à tuberculose dans le groupe de contrôle
  • cneg : le nombre de personnes négatives à tuberculose dans le groupe de contrôle

Nous n’avons pas besoin des autres variables pour le moment.

Nous commençons par calculer les tailles d’effet et leur variance avec la fonction escalc. Pour cela, on précise : - measure= “OR” : pour indiquer que nous souhaitons le rapport de cotes - ai= “tpos” : pour indiquer le nombre de personnes traitées positives - bi= “tneg” : pour indiquer le nombre de personnes traitées négatives - ci= “cpos” : pour indiquer le nombre de contrôles positifs - di= “cneg” : pour indiquer le nombre de contrôles négatifs

Il est à noter qu’il aurait été possible d’utiliser la méthode de Peto (Yusuf et al., 1985) pour cette situation en préciant l’argument measure= “PETO”. Cette méthode est adaptées pour les rapports de cotes. Cette méthode est adaptée pour les intervention ayant de petite tailles d’effets et les événements plutôt rares. En revanche, dans les autres situations, cette méthode est biaisée, c’est pourquoi cette méthode ne peut être recommandée comme approche par défaut dans les méta-analyses. Un des avantages de cette méthode est qu’il n’est pas nécessaire d’appliquer une correction dans les situations où des cellules valent 0 (Higgins et al., 2019).

dat <- escalc(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

Nous observons que le logarithme des rapports de cotes ont été ajoutés au jeu de données.

Données de l'efficacité du vaccin BCG contre la tuberculose après avoir inclu les tailles d'effet

trial

author

year

tpos

tneg

cpos

cneg

ablat

alloc

yi

vi

1

Aronson

1,948

4

119

11

128

44

random

-0.93869414

0.357124952

2

Ferguson & Simes

1,949

6

300

29

274

55

random

-1.66619073

0.208132394

3

Rosenthal et al

1,960

3

228

11

209

42

random

-1.38629436

0.433413078

4

Hart & Sutherland

1,977

62

13,536

248

12,619

52

random

-1.45644355

0.020314413

5

Frimodt-Moller et al

1,973

33

5,036

47

5,761

13

alternate

-0.21914109

0.051951777

6

Stein & Aronson

1,953

180

1,361

372

1,079

44

alternate

-0.95812204

0.009905266

7

Vandiviere et al

1,973

8

2,537

10

619

19

random

-1.63377584

0.227009675

8

TPT Madras

1,980

505

87,886

499

87,892

13

random

0.01202060

0.004006962

9

Coetzee & Berjak

1,968

29

7,470

45

7,232

27

random

-0.47174604

0.056977124

10

Rosenthal et al

1,961

17

1,699

65

1,600

42

systematic

-1.40121014

0.075421726

11

Comstock et al

1,974

186

50,448

141

27,197

18

systematic

-0.34084965

0.012525134

12

Comstock & Webster

1,969

5

2,493

3

2,338

33

systematic

0.44663468

0.534162172

13

Comstock et al

1,976

27

16,886

29

17,825

33

systematic

-0.01734187

0.071635117

Nous pouvons à présent réaliser un modèle à effet aléatoire en utilisant l’algorithme Paule-Mandel en accord avec les recommandatons de Veroniki et al. (2015).

rma.out <- rma(yi = yi, # la colonne dans le jeu de données où on trouve la taille d'effet 
                vi = vi, # la colonne où on trouve l'erreur de mesure  
                data=dat, # le nom du jeu de données 
                method="PM", # utiliser la méthode de Paule-Mandel
                digits=2) # le nombre de décimale à imprimer dans la sortie de résultats
summary(rma.out)
## 
## Random-Effects Model (k = 13; tau^2 estimator: PM)
## 
##   logLik  deviance       AIC       BIC      AICc   
##   -13.10     37.31     30.20     31.33     31.40   
## 
## tau^2 (estimated amount of total heterogeneity): 0.34 (SE = 0.19)
## tau (square root of estimated tau^2 value):      0.58
## I^2 (total heterogeneity / total variability):   92.15%
## H^2 (total variability / sampling variability):  12.73
## 
## Test for Heterogeneity:
## Q(df = 12) = 163.16, p-val < .01
## 
## Model Results:
## 
## estimate    se   zval  pval  ci.lb  ci.ub      
##    -0.75  0.19  -3.99  <.01  -1.11  -0.38  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rappelez-vous que pour les OR, les valeurs calculées avaient subi une transformation logarithmique. Dès lors pour pouvoir un rapport de cote interprétable, il faut appliquer un exponentiel sur la valeur moyenne obtenue dans la sortie de la méta-analyse. On obtient cette valeur, ainsi que l’intervalle de confiance, grâce à la fonction predict. Deux arguments doivent être précisés :

  • object : le modèle de la méta-analyse
  • transf : la transformation à devoir opérer. En l’occurence nous souhaitons un exponentiel. on va préciser “exp”.

Cela nous fournit le rapport de cote moyen obtenu, l’intervalle de confiance à 95% pour la vraie valeur du rapport de cotes moyen et l’intervalle de prédiction.

Encart 1. L’intervalle de confiance et l’intervalle de prédiction.

L’intervalle de confiance correspond aux limites entre lesquelles la taille d’effet moyenne réelle est censée se situer. Concrètement, nous avons 95% de chances que la vraie valeur du rapport de cotes se situe entre ces deux limites.

L’intervalle de prédiction représente l’intervalle à 95% dans lequel nous avons 95% de chances qu’une nouvelle observation se site. L’intervalle de prédiction est généralement plus large que l’intervalle de confiance.

predict(object = rma.out, transf=exp)  
## 
##  pred ci.lb ci.ub pi.lb pi.ub 
##  0.47  0.33  0.68  0.14  1.58

Il ne reste plus qu’à faire le forest plot (voir plus haut pour une description détaillée).

Exercice 2

Réalisez le forest plot en vous appuyant sur celui réalisé dans la partie consacrée aux différences de moyennes standardisées.

Cliquez pour la solution 
par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(rma.out)), "%") # weights in % (rounded)

forest(rma.out, xlim=c(-4,4.5), header=TRUE, top=2, mlab="Summary",
       ilab=dat$weights, ilab.xpos=2.6, ilab.pos=2.6)
text(2.6, -1, "100%", pos=2)
text(2.6, 15 , "Weight", pos=2, font=2)

Exercice 3

Le Syndrome d’Apnées Obstructives du Sommeil (SAOS) occasionne, par ses mécanismes de pression thoracique et de désaturation en oxygène, de plus grands efforts de la part du cœur pour maintenir les fonctions vitales lors des phases du sommeil. Cette charge mécanique entraine une plus grande fragilité de cet organe. À long terme, on peut s’attendre à une augmentation significative des risques cardiovasculaires chez les personnes souffrant de SAOS.

Des auteurs ont souhaité mesurer l’ampleur de ce phénomène en évaluant la survenue d’Accidents Vasculaires Cérébraux (AVC) sur 10 ans après diagnostic de SAOS, ce en comparaison à un échantillon contrôle. Note* : les données présentées dans cet exemple sont fictives

Les données sont présentées dans la feuille “ex2” du fichier excel fourni précédemment.

Réalisez la méta-analyse sur les rapports de cotes et faites le forest plot. Utilisez l’algorithme adapté.

Cliquez pour la solution 
dat <- read_excel("Exercices.xlsx", sheet="ex2")

Les données sont présentées ci-dessous.

Données SAOS - rapports de cotes

study

AVC_SAOS

Safe_SAOS

n_SAOS

AVC_GC

Safe_GC

n_GC

Bae et al., 2021

101

86

187

47

153

200

Heilbrunn et al., 2021

15

6

21

2

15

17

Mashagi et al., 2020

35

37

72

35

50

85

Ottaviani et al., 2020

24

26

50

16

34

50

Xie et al., 2017

12

8

20

5

14

19

Yaggi et al., 2022

17

12

29

8

28

36

##### > Étape n° 4 : Calcul des variances et des tailles d'effet 
dat <- escalc("OR", ai=AVC_SAOS, n1i=n_SAOS, ci=AVC_GC, n2i=n_GC, slab=study, data=dat)
dat
## 
##                    study AVC_SAOS Safe_SAOS n_SAOS AVC_GC Safe_GC n_GC     yi 
## 1       Bae et al., 2021      101        86    187     47     153  200 1.3411 
## 2 Heilbrunn et al., 2021       15         6     21      2      15   17 2.9312 
## 3   Mashagi et al., 2020       35        37     72     35      50   85 0.3011 
## 4 Ottaviani et al., 2020       24        26     50     16      34   50 0.6737 
## 5       Xie et al., 2017       12         8     20      5      14   19 1.4351 
## 6     Yaggi et al., 2022       17        12     29      8      28   36 1.6011 
##       vi 
## 1 0.0493 
## 2 0.8000 
## 3 0.1042 
## 4 0.1720 
## 5 0.4798 
## 6 0.3029
# Modèle à effet fixe (Equal Effect model) 
model_RE <- rma(yi, vi, data=dat, method="PM", digits=2)
model_RE
## 
## Random-Effects Model (k = 6; tau^2 estimator: PM)
## 
## tau^2 (estimated amount of total heterogeneity): 0.39 (SE = 0.40)
## tau (square root of estimated tau^2 value):      0.62
## I^2 (total heterogeneity / total variability):   69.40%
## H^2 (total variability / sampling variability):  3.27
## 
## Test for Heterogeneity:
## Q(df = 5) = 13.61, p-val = 0.02
## 
## Model Results:
## 
## estimate    se  zval  pval  ci.lb  ci.ub      
##     1.19  0.32  3.66  <.01   0.55   1.82  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# transformation des données catégorielles via la fonction exponentielle
# pour s'ajuster au modèle de régression logistique
# À utiliser ++ dans un EE model
predict(model_RE, transf=exp)  
## 
##  pred ci.lb ci.ub pi.lb pi.ub 
##  3.27  1.74  6.17  0.83 12.95
par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(model_RE)), "%") # weights in % (rounded)

forest(model_RE, xlim=c(-3.5,5), header=TRUE, top=2, mlab="Summary",
       atransf=exp, at=log(c(0.125, 0.25, 0.5, 1, 2, 4)), digits=c(2L,3L),
       ilab=dat$weights, ilab.xpos=2.5, ilab.pos=2)
text(2.5, -1, "100%", pos=2)
text(2.5,  8, "Weight", pos=2, font=2)

Les méta-analyses sur des corrélations

Des chercheurs (McDaniel et al., 1994) se sont intéressés à la pertinence des entretiens d’embauche. Plus spécifiquement, ils se sont demandés si la qualité de l’évaluation lors de l’entretien d’embauche était un bon prédicteur de l’adéquation aux tâches professionnelles réalisées par la personne.

Pour obtenir les données, vous pouvez tapez

data(dat.mcdaniel1994)

Les 20 premières lignes du jeu de données sont présentées ci-dessous.

Données de McDaniel et al. (1994) sur le lien entre la qualité de l'entretien d'embauche et l'efficacité sur le poste de travail occupé.

study

ni

ri

type

struct

1

123

0.00

j

s

2

95

0.06

p

u

3

69

0.36

j

s

4

1,832

0.15

j

s

5

78

0.14

j

s

6

329

0.06

j

s

7

153

0.09

j

s

8

29

0.40

j

s

9

29

0.39

s

s

10

157

0.14

s

s

11

149

0.36

s

s

12

92

0.28

j

u

13

15

0.62

j

s

14

15

0.07

j

u

15

170

0.18

j

u

16

19

0.42

j

s

17

19

0.08

j

u

18

68

0.18

p

u

19

93

0.43

j

u

20

57

0.04

j

u

Dans ce jeu de données, les variables sont :

  • study : le numéro de l’étude
  • ni : la taille de l’échantillon
  • ri : la corrélation observée
  • type : le type d’entretien (lié au job, situationnel ou psychologique)
  • struct : la structure de l’entretien (structuré vs. libre)

Il est possible d’obtenir la distribution autour de la corrélation en utilisant la fonction escalc. Il faut dans ce cas préciser que l’argument measure vaut “COR”.

dat <- escalc(measure="COR", ri=ri, ni=ni, slab=study, data=dat.mcdaniel1994)
head(dat)
## 
##   study   ni   ri type struct     yi     vi 
## 1     1  123 0.00    j      s 0.0000 0.0082 
## 2     2   95 0.06    p      u 0.0600 0.0106 
## 3     3   69 0.36    j      s 0.3600 0.0111 
## 4     4 1832 0.15    j      s 0.1500 0.0005 
## 5     5   78 0.14    j      s 0.1400 0.0125 
## 6     6  329 0.06    j      s 0.0600 0.0030

Sachant que les corrélations sont biaisées, il est préférable d’utiliser une transformation z de Fisher avant de réaliser la méta-analyse Fisher & Filon (1898). Appliquez la transformation de Fisher sur les données. Pour cela, l’argument measure doit être “ZCOR”.

Exercice 4

Réalisez la transformation de z de Fisher sur les corrélations des données de McDaniel et al. (1994).Comparez les valeurs des corrélations avec celles des transformations de Fisher.

dat <- escalc(measure="ZCOR", ri=ri, ni=ni, slab=study, data=dat, var.names=c("z.yi", "z.vi"))
head(dat)
## 
##   study   ni   ri type struct     yi     vi   z.yi   z.vi 
## 1     1  123 0.00    j      s 0.0000 0.0082 0.0000 0.0083 
## 2     2   95 0.06    p      u 0.0600 0.0106 0.0601 0.0109 
## 3     3   69 0.36    j      s 0.3600 0.0111 0.3769 0.0152 
## 4     4 1832 0.15    j      s 0.1500 0.0005 0.1511 0.0005 
## 5     5   78 0.14    j      s 0.1400 0.0125 0.1409 0.0133 
## 6     6  329 0.06    j      s 0.0600 0.0030 0.0601 0.0031

L’étape suivante consiste à faire la méta-analyse. Il est techniquement possible de faire la méta-analyse sur les corrélations ou sur les tansformations z de Fisher. Étant donné que la transformation z est moins biaisé, il est préférable d’utiliser la transformation z et d’obtenir le coefficient de corrélation moyen après analyse. Nous avons vu que la méthode la plus adaptée pour des données quantitatives était le maximum de vraisemblance restreint. Nous allons donc utiliser cet algorithme.

model_RE <- rma(yi=z.yi, vi=z.vi, data=dat, method="REML", digits=2)
summary(model_RE)
## 
## Random-Effects Model (k = 160; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc   
##     5.75    -11.50     -7.50     -1.36     -7.42   
## 
## tau^2 (estimated amount of total heterogeneity): 0.03 (SE = 0.00)
## tau (square root of estimated tau^2 value):      0.17
## I^2 (total heterogeneity / total variability):   81.29%
## H^2 (total variability / sampling variability):  5.35
## 
## Test for Heterogeneity:
## Q(df = 159) = 789.73, p-val < .01
## 
## Model Results:
## 
## estimate    se   zval  pval  ci.lb  ci.ub      
##     0.24  0.02  14.00  <.01   0.20   0.27  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La valeur 0.24 est le z de Fisher mais nous souhaiterions plutôt la valeur de la corrélation. Nous pouvons obtenir cette valeur grâce à la fonction predict. Il faut préciser le type de transformation souhaitée dans l’argument transf. En l’occurrence, il s’agit de transformer le z en r ‘transf.ztor’

predict(model_RE, transf=transf.ztor)
## 
##  pred ci.lb ci.ub pi.lb pi.ub 
##  0.23  0.20  0.26 -0.10  0.52

La suite consiste à réaliser le forest plot comme précédemment. La difficulté ici est que le nombre d’études est particulièrement important (càd. 160).

Encart 2. Réaliser un forest plot avec un grand nombre d’études.

Quand la méta-analyse inclut un grand nombre d’études, il existe plusieurs difficultés pour réaliser le forest plot.

  1. la fonction forest renverra un graphique difficilement lisible si vous essayez de le faire directement dans R. Il faut donc l’exporter, par exemple en png.

  2. Il faut trouver les bons ajustements de tailles, à la fois la haut et la largeur. Pour la hauteur, il faut compter au moins 100 pixels par ligne et ajouter une petite marge de manoeuvre. En l’occurrence, nous avons décidé d’opter pour 25000 pour la hauteur. Pour la largeur, il faut qu’elle ne soit pas trop petite pour ne pas que cela paraisse bizarre par rapport à la hauteur. Avoir une largeur de 40% de la taille de la hauteur est un minimum. Nous optons pour 10 000.

  3. Si on réalise le forest plot avec la paramètres ci-dessus sans changer les autres paramètres, on constate qu’il y a superposition entre les intervalles de confiance et le texte. Il faut donc ajuster la position du texte.

png(file='forestplot.png',width= 10000, height=25000, res = 400 ) # Open PNG device with specific file name
par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(model_RE)), "%") 

forest(model_RE, xlim=c(-3.5,5), header=TRUE, top=2, mlab="Summary",  digits=c(2L,3L),
       ilab=dat$weights, ilab.xpos=4.2, ilab.pos=2)
text(4.2, -1, "100%", pos=2)
text(4.2,  162, "Weight", pos=2, font=2)
dev.off()

Exercice 5

D’autres chercheurs se sont intéressés à la relation qui existe entre le nombre d’apnées du sommeil durant la nuit et la qualité de vie des personnes. Des corrélations ont pu être établies entre ces deux variables, mais la littérature demeure encore jeune, avec un risque d’inconsistance majeur : une méta-analyse pourrait permettre d’y voir plus clair dans ces premières données.

Note : les données présentées dans cet exemple sont fictives

Les données sont présentées dans le fichier “Exercice.xlsx” sur la feuille “ex3”. Réalisez un modèle à effet fixe, obtenez la corrélation moyenne, faites le forest plot. Refaites cet exercice en réalisant un modèle aléatoire qui utilise la méthode DerSimonian-Laird.

Cliquez pour obtenir la solution 
dat <- read_excel("Exercices.xlsx", sheet="ex3")

Les données sont présentées ci-dessous.

Données SAOS - Corrélations

study

corr_dep

n

Gottlieb et al., 2019

0.49

4,422

Young et al., 2017

0.60

1,522

Punjabi et al., 2018

0.63

6,294

Yaggi et al., 2005

0.69

1,022

Leao et al., 2016

0.52

73

# Ajout des transformation z de Fisher pour les corrélations et la variance associée
dat <- escalc("ZCOR", ri=corr_dep, ni=n, slab=study, data=dat)
dat
## 
##                   study corr_dep    n     yi     vi 
## 1 Gottlieb et al., 2019     0.49 4422 0.5361 0.0002 
## 2    Young et al., 2017     0.60 1522 0.6931 0.0007 
## 3  Punjabi et al., 2018     0.63 6294 0.7414 0.0002 
## 4    Yaggi et al., 2005     0.69 1022 0.8480 0.0010 
## 5     Leao et al., 2016     0.52   73 0.5763 0.0143
# Réalisation d'un modèle à effets fixes
model_EE <- rma(yi, vi, data=dat, method="EE", digits=2)
model_EE
## 
## Equal-Effects Model (k = 5)
## 
## I^2 (total heterogeneity / total variability):   97.24%
## H^2 (total variability / sampling variability):  36.18
## 
## Test for Heterogeneity:
## Q(df = 4) = 144.72, p-val < .01
## 
## Model Results:
## 
## estimate    se   zval  pval  ci.lb  ci.ub      
##     0.68  0.01  77.90  <.01   0.66   0.69  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Transformation du z de Fisher en coefficient de corrélation 
# calcul de l'intervalle de confiance et de l'intervalle de prédiction 
predict(model_EE, transf=transf.ztor)
## 
##  pred ci.lb ci.ub 
##  0.59  0.58  0.60
# Réalisation du forest plot

par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(model_EE)), "%") # weights in % (rounded)

forest(model_EE, xlim=c(-2,4), header=TRUE, top=2, mlab="Summary",
       atransf=transf.ztor, at=transf.rtoz(c(-0.5, 0, 0.5, 0.85)),
       ilab=dat$weights, ilab.xpos=2.2, ilab.pos=2)
text(2.2, -1, "100%", pos=2)
text(2.2,  7, "Weight", pos=2, font=2)

# Modèle aléatorisé (Random Effect model) avec la méthode DerSimonian-Laird
model_RE <- rma(yi, vi, data=dat, method="DL", digits=2)
model_RE
## 
## Random-Effects Model (k = 5; tau^2 estimator: DL)
## 
## tau^2 (estimated amount of total heterogeneity): 0.02 (SE = 0.02)
## tau (square root of estimated tau^2 value):      0.13
## I^2 (total heterogeneity / total variability):   97.24%
## H^2 (total variability / sampling variability):  36.18
## 
## Test for Heterogeneity:
## Q(df = 4) = 144.72, p-val < .01
## 
## Model Results:
## 
## estimate    se   zval  pval  ci.lb  ci.ub      
##     0.69  0.06  11.32  <.01   0.57   0.81  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict(model_RE, transf=transf.ztor)
## 
##  pred ci.lb ci.ub pi.lb pi.ub 
##  0.60  0.51  0.67  0.39  0.75
par(mar=c(4,4,2,2))

dat$weights <- paste0(round(weights(model_RE)), "%") # weights in % (rounded)

forest(model_RE, xlim=c(-2,4), header=TRUE, top=2, mlab="Summary",
       atransf=transf.ztor, at=transf.rtoz(c(-0.5, 0, 0.5, 0.85)),
       ilab=dat$weights, ilab.xpos=2.2, ilab.pos=2)
text(2.2, -1, "100%", pos=2)
text(2.2,  7, "Weight", pos=2, font=2)

Codes pour les packages meta et dmetar à venir

Références

Alexander, R. A., Scozzaro, M. J., & Borodkin, L. J. (1989). Statistical and empirical examination of the chi-square test for homogeneity of correlations in meta-analysis. Psychological Bulletin, 106(2), 329–331. https://doi.org/10.1037/0033-2909.106.2.329
Balduzzi, S., Rücker, G., & Schwarzer, G. (2019). How to perform a meta-analysis with r: A practical tutorial. Evidence Based Mental Health, 22(4), 153–160. https://doi.org/10.1136/ebmental-2019-300117
Ben-Shachar, M. S., Lüdecke, D., & Makowski, D. (2020). effectsize: Estimation of effect size indices and standardized parameters. Journal of Open Source Software, 5(56), 2815. https://doi.org/10.21105/joss.02815
Borenstein, M., Hedges, L. V., Higgins, J. P. T., & Rothstein, H. R. (2009). Introduction to meta‐analysis. Wiley. https://doi.org/10.1002/9780470743386
Fisher, K., & Filon, L. N. G. (1898). Proceedings of the Royal Society of London, 62(379–387), 173–176. https://doi.org/10.1098/rspl.1897.0091
Furukawa, T. A., McGuire, H., & Barbui, C. (2003). Low dosage tricyclic antidepressants for depression. Cochrane Database of Systematic Reviews. https://doi.org/10.1002/14651858.cd003197
Harrer, M., Cuijpers, P., A, F. T., & Ebert, D. D. (2021). Doing meta-analysis with R: A hands-on guide (1st ed.). Chapman & Hall/CRC Press.
Harrer, M., Cuijpers, P., Furukawa, T., & Ebert, D. D. (2019). Dmetar: Companion r package for the guide ’doing meta-analysis in r’. http://dmetar.protectlab.org/
Hedges, L., & Olkin, I. (1985). Statistical methods for meta-analysis. Elsevier. https://doi.org/10.1016/c2009-0-03396-0
Higgins, J. P. T., Li, T., & Deeks, J. (2019). Choosing effect measures and computing estimates of effect. In J. P. T. Higgins, J. Thomas, J. Chandler, M. Cumpston, T. Li, M. J. Page, & V. A. Welch (Eds.), Cochrane handbook for systematic reviews of interventions. (pp. 143–176). Wiley.
Laird, N. M., & Mosteller, F. (1990). Some statistical methods for combining experimental results. International Journal of Technology Assessment in Health Care, 6(1), 5–30. https://doi.org/10.1017/s0266462300008916
McDaniel, M. A., Whetzel, D. L., Schmidt, F. L., & Maurer, S. D. (1994). The validity of employment interviews: A comprehensive review and meta-analysis. J. Appl. Psychol., 79(4), 599–616.
Poole, C., & Greenland, S. (1999). Random-effects meta-analyses are not always conservative. American Journal of Epidemiology, 150(5), 469–475. https://doi.org/10.1093/oxfordjournals.aje.a010035
Schwarzer, G., Carpenter, J. R., & Rücker, G. (2015). Meta-analysis with r. Springer International Publishing. https://books.google.fr/books?id=HQO0CgAAQBAJ
Stanley, T. D., Doucouliagos, H., & Ioannidis, J. P. A. (2022). Beyond random effects: When small-study findings are more heterogeneous. Advances in Methods and Practices in Psychological Science, 5(4), 251524592211204. https://doi.org/10.1177/25152459221120427
Veroniki, A. A., Jackson, D., Viechtbauer, W., Bender, R., Bowden, J., Knapp, G., Kuss, O., Higgins, J. P., Langan, D., & Salanti, G. (2015). Methods to estimate the between‐study variance and its uncertainty in meta‐analysis. Research Synthesis Methods, 7(1), 55–79. https://doi.org/10.1002/jrsm.1164
Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03
Yusuf, S., Peto, R., Lewis, J., Collins, R., & Sleight, P. (1985). Beta blockade during and after myocardial infarction: An overview of the randomized trials. Progress in Cardiovascular Diseases, 27(5), 335–371. https://doi.org/10.1016/s0033-0620(85)80003-7