Auteur : Nicolas Stefaniak

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

Prérequis

Avoir installé et chargé les packages metafor (Viechtbauer & Cheung, 2010), meta (Balduzzi et al., 2019) et dmetar (Harrer et al., 2019a).

Si les packages sont installés, il suffit d’utiliser la fonction require avec le nom du package.

require(metafor)
require(dmetar)
require(meta)

Si ‘metafor’ou ’mta’ n’est/ne sont pas installé(s), avant d’utiliser la fonction require, il faut taper :

install.packages("metafor")
install.packages("meta")

Si ‘dmetar’ n’est pas installé, pour l’installer, il faut taper :

if (!require("remotes")) {
  install.packages("remotes")
}
remotes::install_github("MathiasHarrer/dmetar")

introduction

Pour savoir s’il est pertinent ou non d’aller explorer la présence possible de valeurs influentes, il n’existe pas de règles absolues, mais plutôt des heuristiques qui, comme toutes les heuristiques ont leur propres limites et sont arbitraires. Dans le cas des méta-analyses, l’heuristique consiste à aller explorer la possible présence de valeurs influentes quand la valeur de \(I^2\) est supérieure à 50%. Pour Harrer (2021), cette règle arbitraire a au moins l’avantage de pouvoir préciser a priori les conditions qui vont amener à aller explorer la présence de valeurs influentes, et de pouvoir le faire de manière consistante.

Le point critique est de supprimer des valeurs en considérant qu’elles sont influentes car les résultats de l’étude ne nous plaisent pas. Il faut toujours avoir une sérieuse raison pour le faire. Le danger est qu’on puisse le faire de manière inconsciente.

En revanche, comme l’estimation des résultats peut être considérablement biaisée à cause d’une ou de quelques études qui ont des tailles d’effet extrêmes, cela peut être une bonne idée de relancer l’analyse une fois que ces valeurs influentes ont été supprimées. Il s’agit aussi de déterminer si les résultats ne dépendent pas de manière trop conséquente d’une seule étude.

Si nous observons la Figure 1, il apparaît sur le graphique qu’une étude (l’étude 7) semble se distinguer de manière considérable de la manière dont les autres études se comportent. La question qui se pose est de savoir comment il est possible d’objectiver que cette étude se comporte de manière considérablement différente. Pour cela, il existe plusieurs méthodes (Viechtbauer & Cheung, 2010).

Figure 1. Forest plot pédagogique avec une valeur influente

une approche basique

Une première manière pour identifier les valeurs influentes est de considérer qu’une étude est influente si les limites de son intervalle de confiance ne tombent pas dans les limites de l’intervalle de confiance de la moyenne pondérée des tailles d’effet. En d’autres termes, la taille d’effet est une valeur aberrante si sa valeur est tellement extrême qu’elle diffère significativement de la taille d’effet global.

Dans l’exemple de la Figure 1, il y a 7 études fictives pour lesquelles nous avons décidé d’avoir la même erreur de mesure , ce qui impliquerait u,e taille d’échantillon strictement identique pour chaque étude. Ce choix a été fait pour faciliter les calculs. Les 7 études présentées dans la Figure 1 reportent les tailles d’effet suivantes :

## [1]  0.40  0.50  0.60  0.66  0.62  0.52 -1.20

La moyenne de ces tailles d’effet vaut :

mean(mean1)
## [1] 0.3

Sachant que nous avons fixé l’écart-type autour de chaque taille d’effet à 0.42, on peut calculer la pondération de chaque étude, la variance globale et donc l’écart-type global de la distribution des tailles d’effet1 :

# la pondération d'une étude 
1/(0.42^2) 
## [1] 5.668934
# La variance vaut donc 1 divisé par 7 fois cette pondération
var<-1/(7*(1/(0.42^2)) )
var
## [1] 0.0252
# et en faisant la racine carré nous obtenons l'écart-type
sqrt(var)
## [1] 0.1587451

Il s’ensuit que l’intervalle de confiance vaut :

0.3-1.96*sqrt(var) # limite inférieure 
## [1] -0.01114035
0.3+1.96*sqrt(var) # limite supérieure
## [1] 0.6111404

Si l’intervalle de confiance d’une étude ne tombe pas à l’intérieur de ces limites, alors c’est un outlier.

En l’occurrence, l’intervalle de confiance autour de l’étude dpnt la taille d’effet est à -1.2 (l’étude 7) est :

-1.2-1.96*0.42 # limite inférieure 
## [1] -2.0232
-1.2+1.96*.42 # limite supérieure
## [1] -0.3768

Il s’agit d’une valeur influente qu’il faudrait supprimer comme on avait pu l’observer assez aisément sur le forest plot voir Figure 1.

Cette méthode est relativement peu utilisée et nous verrons plus loin des méthodes plus efficaces et robustes pour détecter les valeurs influentes/aberrantes. Néanmoins, elle existe et il peut être utile d’illustrer comment l’utiliser de manière concrète. Le package ‘metafor’ n’a pas de fonction qui utilise cette méthode, mais elle est préente dans “dmetar” (Harrer et al., 2019a). Nous allons donc utiliser ce package.

Il est à noter que ‘dmetar’ est compatible avec le package ‘meta’ (Balduzzi et al., 2019). Nous allons donc réaliser la méta-analyse avec les deux packages en choisissant les paramètres qui vont permettre de répliquer l’analyse réalisée par Harrer et al. (2021), chapitre 4.

Pour notre exemple de travail, nous pouvons utiliser les données “ThirdWave” issues du package ‘dmetar’. Ces données sont des tailles d’effet sur l’efficacité des thérapies cognitivo-comportementales de troisième génération pour lutter contre le stress perçu chez des étudiants par rapport à un groupe de contrôle inactif.

Les tailles d’effets sont des différences de moyennes standardisées voir Tableau 1.

Données sur l'efficacité des TCC de 3e génération entre un groupe traité et un groupe de contrôle inactif

Author

TE

seTE

RiskOfBias

TypeControlGroup

InterventionDuration

InterventionType

ModeOfDelivery

Call et al.

0.71

0.26

high

WLC

short

mindfulness

group

Cavanagh et al.

0.35

0.20

low

WLC

short

mindfulness

online

DanitzOrsillo

1.79

0.35

high

WLC

short

ACT

group

de Vibe et al.

0.18

0.12

low

no intervention

short

mindfulness

group

Frazier et al.

0.42

0.14

low

information only

short

PCI

online

Frogeli et al.

0.63

0.20

low

no intervention

short

ACT

group

Gallego et al.

0.72

0.22

high

no intervention

long

mindfulness

group

Hazlett-Stevens & Oren

0.53

0.21

low

no intervention

long

mindfulness

book

Hintz et al.

0.28

0.17

low

information only

short

PCI

online

Kang et al.

1.28

0.34

low

no intervention

long

mindfulness

group

Kuhlmann et al.

0.10

0.19

high

no intervention

short

mindfulness

group

Lever Taylor et al.

0.39

0.23

low

WLC

long

mindfulness

book

Phang et al.

0.54

0.24

low

no intervention

short

mindfulness

group

Rasanen et al.

0.43

0.26

low

WLC

short

ACT

online

Ratanasiripong

0.52

0.35

high

no intervention

short

mindfulness

group

Shapiro et al.

1.48

0.32

high

WLC

long

mindfulness

group

Song & Lindquist

0.61

0.23

high

WLC

long

mindfulness

group

Warnecke et al.

0.60

0.25

low

information only

long

mindfulness

online

Il est à noter que dans ce jeu de données, la mesure de dispersion est déjà présente. Il ’agit de l’erreur-type et non de la variance.

Nous commencerons par réaliser la méta-analyse avec le package ‘meta’. Pour cela, il faut utiliser la fonction metagen :

  • TE : les tailles d’effets

  • seTE : la mesure de dispersion sur les tailles d’effet

  • studlab : la colonne qui contient le nom des études

  • data : le nom du jeu de données

  • sm : le type de mesure standardisée. Ici, il s’agit d’une différence de moyenne standardisée. Les autres possibles sont “RD” (différence de risques), “RR” (risque relatif), “OR” (odd ratio), “HR” (ratio de hasard), “ASD” (la transformation arcinus des probabilités), “MD” (la différence de moyennes), et “ROM (le rapport des moyennes)

  • fixed : logique. Faut-il calculer un modèle à effets fixes ?

  • random : logique. Faut-il calculer un modèle à effets aléatoires ?

  • method.tau : la méthode pour calculer la variance aléatoire. Nous utiliserons ici le maximum de vraisemblance restreint.

  • method.random.ci : la méthode pour calculer les tests et l’intervalle de confiance de l’effet aléatoire. En l’occurrence, c’est la méthode Hartung-Knapp (Knapp & Hartung, 2003).

  • title : le titre de l’analyse, qui sera utile non seulement pour les résultats mais aussi pour le forest plot.

m.gen <- metagen(TE = TE,
                 seTE = seTE,
                 studlab = Author,
                 data = ThirdWave,
                 sm = "SMD",
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 method.random.ci = "HK",
                 title = "Third Wave Psychotherapies")
summary(m.gen)
## Review:     Third Wave Psychotherapies
## 
##                           SMD            95%-CI %W(random)
## Call et al.            0.7091 [ 0.1979; 1.2203]        5.0
## Cavanagh et al.        0.3549 [-0.0300; 0.7397]        6.3
## DanitzOrsillo          1.7912 [ 1.1139; 2.4685]        3.8
## de Vibe et al.         0.1825 [-0.0484; 0.4133]        7.9
## Frazier et al.         0.4219 [ 0.1380; 0.7057]        7.3
## Frogeli et al.         0.6300 [ 0.2458; 1.0142]        6.3
## Gallego et al.         0.7249 [ 0.2846; 1.1652]        5.7
## Hazlett-Stevens & Oren 0.5287 [ 0.1162; 0.9412]        6.0
## Hintz et al.           0.2840 [-0.0453; 0.6133]        6.9
## Kang et al.            1.2751 [ 0.6142; 1.9360]        3.9
## Kuhlmann et al.        0.1036 [-0.2781; 0.4853]        6.3
## Lever Taylor et al.    0.3884 [-0.0639; 0.8407]        5.6
## Phang et al.           0.5407 [ 0.0619; 1.0196]        5.3
## Rasanen et al.         0.4262 [-0.0794; 0.9317]        5.1
## Ratanasiripong         0.5154 [-0.1731; 1.2039]        3.7
## Shapiro et al.         1.4797 [ 0.8618; 2.0977]        4.2
## Song & Lindquist       0.6126 [ 0.1683; 1.0569]        5.7
## Warnecke et al.        0.6000 [ 0.1120; 1.0880]        5.2
## 
## Number of studies: k = 18
## 
##                              SMD           95%-CI    t  p-value
## Random effects model (HK) 0.5771 [0.3782; 0.7760] 6.12 < 0.0001
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0820 [0.0295; 0.3533]; tau = 0.2863 [0.1717; 0.5944]
##  I^2 = 62.6% [37.9%; 77.5%]; H = 1.64 [1.27; 2.11]
## 
## Test of heterogeneity:
##      Q d.f. p-value
##  45.50   17  0.0002
## 
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Hartung-Knapp adjustment for random effects model (df = 17)

Nous allons également faire cette analyse avec le package ‘metafor’ (Viechtbauer & Cheung, 2010). Afin de reproduire les résultats de Harrer et al. (2021), chapitre 4, nous utiliserons les mêmes paramètres que ces derniers avec le package ‘metafor’. Ainsi, nous fixerons l’argument method=“REML”. Pour préciser la manière de calculer l’intervalle de confiance et les tests, il faut utiliser l’argument test qui prendra la valeur knha pour s’assurer que ce soit la méthod Hartung-Knapp (Knapp & Hartung, 2003). Pour rappel, ‘metafor’ requiert normalement que la mesure de dispersion des tailles d’effet soit les variances autour des tailles d’effet. Cependant, dans notre jeu de données, nous ne disposons pas de cette information mais nous avons les erreurs-types. Cela ne présente pas de difficulté pour la fonction rma qui a un argument sei qu’on peut utiliser à la place de l’argument vi. Pour la suite du chapitre, nous allons également préciser la colonne où se situe le nom des études avec l’argument slab.

rma.out<-rma(yi = TE, sei = seTE, data =ThirdWave, method="REML" ,test="knha", slab=Author)

rma.out
## 
## Random-Effects Model (k = 18; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0820 (SE = 0.0458)
## tau (square root of estimated tau^2 value):      0.2863
## I^2 (total heterogeneity / total variability):   64.63%
## H^2 (total variability / sampling variability):  2.83
## 
## Test for Heterogeneity:
## Q(df = 17) = 45.5026, p-val = 0.0002
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub      
##   0.5771  0.0943  6.1222  17  <.0001  0.3782  0.7760  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Identifiez que les résultats sont pratiquement exactement identiques à ceux présentés précédemment, ce qui nous permettra de continuer avec ‘metafor’ pour la suite.

À présent que nous avons réalisé la méta-analyse, nous pouvons aller identifier la présence des outliers avec la méthode que nous venons de décrire. Ceci est possible avec la fonction find.outliers du package ‘dmetar’. Cette fonction a comme argument la sortie de résultats de la métaanalyse réalisée avec le package ‘meta’, à savoir l’objet code>metagen.

outlier.out<-find.outliers(m.gen)
outlier.out
## Identified outliers (random-effects model) 
## ------------------------------------------ 
## "DanitzOrsillo", "Shapiro et al." 
##  
## Results with outliers removed 
## ----------------------------- 
## Review:     Third Wave Psychotherapies
## 
## Number of studies: k = 16
## 
##                              SMD           95%-CI    t  p-value
## Random effects model (HK) 0.4528 [0.3257; 0.5800] 7.59 < 0.0001
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0139 [0.0000; 0.1032]; tau = 0.1180 [0.0000; 0.3213]
##  I^2 = 24.8% [0.0%; 58.7%]; H = 1.15 [1.00; 1.56]
## 
## Test of heterogeneity:
##      Q d.f. p-value
##  19.95   15  0.1739
## 
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Hartung-Knapp adjustment for random effects model (df = 15)

Cette fonction identife que deux études peuvent être considées comme des outliers (“DanitzOrsillo” et “Shapiro et al.” ) et relance l’analyse sans ces études. On sait que ces deux études ont été exclues car, si on demande le poids des études, on se rend compte qu’ils sont fixés à 0.

poids<-data.frame(etude = outlier.out$m.random$studlab, outlier.out$m.random$w.random)
myft<-flextable(poids )
  myft<-set_table_properties(myft, layout = "autofit")
  myft<-  style(myft,
        pr_p = fp_par(
          text.align = "center", padding = 1))
  myft <- colformat_double(
  x = myft,
  big.mark = ",", digits = 2, na_str = "N/A"
)
myft<- set_caption(myft, caption = "Poids des études après avoir supprimé les valeurs identifiées comme influentes par la fonction find.outliers")                
myft
Poids des études après avoir supprimé les valeurs identifiées comme influentes par la fonction find.outliers

etude

outlier.out.m.random.w.random

Call et al.

12.20

Cavanagh et al.

19.06

DanitzOrsillo

0.00

de Vibe et al.

35.99

Frazier et al.

28.67

Frogeli et al.

19.11

Gallego et al.

15.53

Hazlett-Stevens & Oren

17.18

Hintz et al.

23.73

Kang et al.

7.84

Kuhlmann et al.

19.29

Lever Taylor et al.

14.89

Phang et al.

13.59

Rasanen et al.

12.43

Ratanasiripong

7.28

Shapiro et al.

0.00

Song & Lindquist

15.31

Warnecke et al.

13.17

Remarquez que la taille d’effet moyenne a diminué et vaut à présent 0.45 au lieu de 0.5571. De même, l’intervalle de confiance est beaucoup plus précis après avoir supprimé les outliers. Pour comprendre ce phénomène, vous pouvez vous reporter au chapitre consacré aux valeurs influentes. De même les mesures d’hétérogénéité ont considérablement diminuées : le \(I^2\) passe de 63% à 25%, l’intervalle de confiance de \(\tau^2\) inclut le 0 et la statistique Q est devenue non significative, ce qui signifie qu’on ne peut plus exclure l’idée qu’il existe une taille d’effet commune à toutes les études.

Les analyses d’influence

À présent que nous avons examiné une méthode assez basique, nous pouvons commencer à examiner des méthodes un peu plus subtiles d’identifier les valeurs susceptibles d’influencer les résultats. En effet, une étude peut avoir une grande influence, même lorsque sa taille d’effet n’est pas particulièrement élevée ou faible (Harrer et al., 2021). C’est notamment le cas, lorsque la taille d’effet globale est considérée comme significative en raison d’une étude avec un échantillon particulièrement grand. Cela signifie donc que les résultats ne sont pas robustes si le seul fait d’enlever une étude de la méta-analyse amène à une conclusion opposée. Il est dès lors nécessaire d’en avoir conscience et de pouvoir en informer les lecteurs.

Harrer et al. (Harrer et al., 2021) soulignent que les outliers n’ont pas nécessairement un impact considérable sur les résultats(mais peuvent en avoir un et c’est souvent le cas) alors que les valeurs influentes ont pour caractéristique de modifier les résultats.

Ainsi, il est nécessaire d’utiliser des outils permettant de diagnostiquer la présence de valeurs influentes (Viechtbauer & Cheung, 2010)

Le graphique de Baujat

Une première manière d’évaluer l’influence d’une étude est de réaliser un graphique de Baujat (Baujat et al., 2002) pour avoir une représentation graphique des études qui pourraient avoir une influence sur les résultats. Ce graphique permet d’identifier de manière visuelle à quel point une étude contribue à l’hétérogénéité. L’axe des abcisses représente la contribution à l’hétérogénéité et l’axe des ordonnées représente l’influence sur la taille d’effet globale. L’influence est calculée en refaisant l’analyse autant de fois qu’il y a d’études en supprimant à chaque fois une étude, selon le principe du ‘leave-one-out’. Ainsi, on peut avoir la différence entre le modèle avec et sans l’étude. Les études les plus à droites sont celles qui contribuent le plus à l’hétérogénéité et celles qui sont les plus en haut sont celles qui contribue le plus à la taille d’effet globale (voir Figure 2).

Dans le package ‘metafor’, le graphique de Baujat est amélioré pour que l’axe des x corresponde aux résidus de Pearson au carré et l’axe des y corresponde à la différence standardisée au carré entre les valeur prédites et ajustée (Viechtbauer, 2020). L’interprétation reste néanmoins sensiblement la même.

baujat(rma.out, symbol="slab")

Figure 2. Graphique de Baujat des mesures d’influences dans une méta-analyse

Obtenir le graphique de Baujat avec le package ‘meta’

On peut également obtenir ce graphique avec le package ‘meta’. Notez qu’il y a de petites différences par rapport au graphique précédent car la procédure n’est pas exactement la même.

meta::baujat(m.gen)

Figure 3. Graphique de Baujat avec le pcakge “meta”

Comme précédemment, ce sont les deux mêmes études qui semblent avoir la plus grande influence sur les résultats.

Les mesures d’influence

Dans le package ‘metafor’, toutes les mesures d’influence peuvent être obtenues par la fonction influence. Il suffit d’indiquer le nom de l’objet dans lequel nous avons stocké les résultats de la méta-analyse obtenus avec la fonction rma (voir Tableau 3).

inf.out<-influence(rma.out)
Mesures d'influences de la méta-analyse des TCC de 3ème génération.

Etude

rstudent

dffits

cook.d

cov.r

tau2.del

QE.del

hat

weight

inf

Call et al.

0.31

0.04

0.00

1.13

0.09

44.71

0.05

5.04

Cavanagh et al.

-0.60

-0.19

0.04

1.12

0.09

45.07

0.06

6.27

DanitzOrsillo

3.01

0.86

0.55

0.63

0.04

30.82

0.04

3.75

*

de Vibe et al.

-1.26

-0.34

0.11

1.04

0.07

37.74

0.08

7.88

Frazier et al.

-0.46

-0.18

0.03

1.15

0.09

45.32

0.07

7.34

Frogeli et al.

0.13

-0.02

0.00

1.15

0.09

44.88

0.06

6.27

Gallego et al.

0.36

0.05

0.00

1.13

0.09

44.26

0.06

5.70

Hazlett-Stevens & Oren

-0.14

-0.09

0.01

1.15

0.09

45.45

0.06

5.98

Hintz et al.

-0.84

-0.25

0.06

1.10

0.09

44.01

0.07

6.85

Kang et al.

1.54

0.38

0.14

0.93

0.07

39.83

0.04

3.86

Kuhlmann et al.

-1.33

-0.31

0.09

1.01

0.07

41.50

0.06

6.30

Lever Taylor et al.

-0.48

-0.16

0.03

1.12

0.09

45.34

0.06

5.59

Phang et al.

-0.10

-0.07

0.01

1.14

0.09

45.44

0.05

5.33

Rasanen et al.

-0.37

-0.13

0.02

1.13

0.09

45.46

0.05

5.09

Ratanasiripong

-0.13

-0.06

0.00

1.11

0.09

45.49

0.04

3.68

Shapiro et al.

2.22

0.65

0.35

0.79

0.05

35.21

0.04

4.16

Song & Lindquist

0.08

-0.03

0.00

1.14

0.09

45.15

0.06

5.66

Warnecke et al.

0.04

-0.03

0.00

1.14

0.09

45.26

0.05

5.25

Remarquez que dans le tableau 3, la colonne “inf” indique par une astérisque les études qui remplissent au moins un des critères pour être considérée comme influente. Il faut noter que, dès lors qu’une observation est considérée comme influente sur une seule mesure, elle sera identifiée comme telle sur l’ensemble des mesures.

Nous allons à présent nous intéresser à chacune de ces mesures et des outils d’aide au diagnostic que nous pouvons utiliser. Il faut savoir qu’il est possible de faire une représentation graphique de l’ensemble des indicateurs avec la fonction plot en précisant qu’il faut représenter graphiquement la sortie de résultats de la fonction influence. Il est possible de customiser le graphique mais cela ne semble pas être le point le plus important ici. Ainsi, la décision de ne pas perdre le lecteur dans des détails techniques a été prise. La Figure 4 représente 8 mesures de valeurs influentes de la méta-analyse.

plot(inf.out)

Figure 4. Huit mesures de diagnostic d’influence

Les résidus studentisés

Il s’agit de l’écart entre la taille d’effet estimée de chaque étude et la taille d’effet globale. Ces résidus sont studentisés.

On obtient cette valeur en calculant l’effet global en ayant enlevé une étude selon le principe de la méthode “en laisser une hors”. Les résidus qui en résultent sont ensuite standardisés par la variance externe à cet effet, le \(\tau^2\) estimé à partir de la taille d’effet combinée externe à cet effet et la variance sans k (Equation 1).

\[t_k =\frac{\hat\theta-\hat\mu_{\text{\k}}}{\sqrt{Var(\mu_\text{\k})+\tau^2_\text{\k}+s^2_\text{k}}}\text{ (1)}\]

En faisant l’hypothèse qu’une étude k donnée s’ajuste correctement à la méta-analyse, les trois termes du dénominateur vont capturer les sources de variabilité qui déterminent à quel point une taille d’effet diffère de la taille d’effet moyenne. Ces sources de variabilité sont :

  1. l’erreur d’échantillonnage de l’étude k ;
  2. la variance des vraies tailles d’effet ;
  3. l’imprécision de la taille d’effet moyenne.

Si une étude ne s’ajuste pas correctement, on peut s’attendre à ce que le résidu soit plus large qu’attendu, amenant à des valeurs de t plus grandes, ce qui indique que l’étude est possiblement influente.

Nous pouvons examiner ces valeurs plus en détails dans la Figure 5.

plot(inf.out, plotinf=1)

Figure 5. Résidus studentisés

Les paramètres par défaut dans le package ‘metafor’ (Viechtbauer & Cheung, 2010) utiliser comme critère une valeur absolue supérieure à 3 pour considérer qu’un résidu studentisé indique la présence d’une valeur influente. Nous voyons ici que l’étude 3, c’est-à-dire celle de DanitzOrsillo est considérée comme une valeur influente.

La valeur DFFITS

La valeur de DFFIT est très similaire à celle du résidu studentisé. On se rend compte que sa forme est sensiblement la même. Cet indicateur montre à quel point les effets combinés changent lorsqu’une étude k est supprimée, exprimée en écart-type. À nouveau, les plus grandes valeurs indiquent que la valeur est potentiellement influente.

Nous pouvons examiner ces valeurs plus en détails dans la Figure 6.

plot(inf.out, plotinf=2)

Figure 6. Graphique des DFFITS.

Une valeur est considérée comme influente si son DFFITS est supérieur à la valeur obtenue par l’équation 2:

\[3 \times \sqrt{\frac{p}{k-p}}\text{ (2)}\] où p représente le nombre de prédicteurs (un seul en l’occurrence) et k est le nombre d’études. Concrètement, cela donne :

\[3 \times \sqrt{\frac{1}{18-1}}=0.72 \] Ainsi, toutes les études dont la valeur de DFFIT est supérieure à 0.72 vont être considérées comme influentes. C’est à nouveau la même étude qui est considérée comme influente.

La distance de Cook

La distance de Cook (\(D_k\)) d’une étude est obtenue par l’équation 3 (Harrer et al., 2021) :

\[ \frac{(\hat\mu-\hat\mu_{\text{\k}})^2}{\sqrt{s^2_k+\tau^2}}\text{ (3)}\]

où \(\hat\mu\) est la moyenne pondérée de toutes les tailles d’effet, et \(\hat\mu_{\text{\k}}\) est la moyenne pondérée des tailles d’effet sans l’étude k. Les autres paramètres de l’équation ont déjà été définis.

Nous pouvons examiner plus en détails les valeurs de la distance de Cook dans la Figure 7.

plot(inf.out, plotinf=3)

Figure 7. Graphique de la distance de Cook

Plus une valeur est grande, plus une observation va être considérée comme influente. Concrètement, une valeur va être considérée comme influente si, pour p degrés de liberté, la distribution \(\chi^2\) est inférieure à 0.50. Dans notre exemple, \(p = 1\). En l’occurrence, le cut-off est donc de

qchisq(.5, # une probabilité de 0.50. 
       1) # le nombre de degrés de liberté 
## [1] 0.4549364

Ainsi, toute distance de Cook supérieure à 0.455 va être considérée comme influente. Dans notre exemple, c’est encore la même étude dont la distance de Cook est supérieure à 0.455.

Le rapport de covariance

Le rapport de covariance d’une étude k est obtenu en divisant la variance de la taille d’effet combinée, à savoir son erreur-type au carré, avec et sans l’étude k (Viechtbauer & Cheung, 2010), qu’on obtient en utilisant l’équation 4.

\[\text{Rapport de Covariance} = \frac{Var(\hat\mu_{\text{\k}})}{Var(\hat\mu)}\text{ (4)}\]

Si le rapport est inférieur à 1, cela signifie que l’estimation de variance est plus précise lorsque l’étude n’est pas présente dans les analyses.

Nous pouvons examiner plus en détails les valeurs du rapport de covariance dans la Figure 8.

plot(inf.out, plotinf=4)

Figure 8. Graphique du rapport de covariance.

Il n’existe cependant pas de règle claire sur la manière d’interpréter cette valeur. Néanmoins, on observe sur la graphique que l’étude 3 est celle pour laquelle la suppression de l’étude amènerait à une réduction substantielle de la variance réiduelle.

La méthode “leave one out” (d’exclusion d’une valeur) sur la valeur de \(T^2\) et de Q

Les Figure 9 et 10 indiquent l’évolution de l’hétérogénéité, estimé par le \(\tau^2\) (Figure 9) et par le Q (Figure 10), lorsque ces valeurs sont calculées avec ou sans l’étude k.

plot(inf.out, plotinf=5)

Figure 9. Changement de l’hétérogénéité, estimée par le \(\tau^2\) par suppression d’une étude.

plot(inf.out, plotinf=6)

Figure 10. Changement de l’hétérogénéité, estimée par la statistique \(Q\) par suppression d’une étude.

Comme il est désirable que l’hétérogénéité soit faible, les diminutions substantielle du Q, et surtout du \(\tau^2\) est souhaité puisque cela indique une plus faible hétérogénéité.

Quelques options spécifiques au package meta et dmetar

En réalité, il est possible d’utiliser cette méthode sur n’importe quel estimateur. Le package meta permet d’obtenir les résultats sur la taille d’effet moyenne estimée et sur le \(I^2\) avec la fonction metainf. Il suffit d’indiquer l’objet issu de la fonction metagen que nous avons utilisée précédemment. En l’occurrence, l’objet s’appelait m.gen.

metainf(m.gen)
## Review:     Third Wave Psychotherapies
## 
## Influential analysis (random effects model (HK))
## 
##                                    SMD           95%-CI  p-value   tau^2
## Omitting Call et al.            0.5736 [0.3616; 0.7855] < 0.0001  0.0902
## Omitting Cavanagh et al.        0.5960 [0.3843; 0.8077] < 0.0001  0.0916
## Omitting DanitzOrsillo          0.5074 [0.3486; 0.6663] < 0.0001  0.0354
## Omitting de Vibe et al.         0.6077 [0.4040; 0.8114] < 0.0001  0.0738
## Omitting Frazier et al.         0.5946 [0.3804; 0.8088] < 0.0001  0.0949
## Omitting Frogeli et al.         0.5787 [0.3644; 0.7931] < 0.0001  0.0941
## Omitting Gallego et al.         0.5718 [0.3593; 0.7843] < 0.0001  0.0904
## Omitting Hazlett-Stevens & Oren 0.5855 [0.3714; 0.7996] < 0.0001  0.0946
## Omitting Hintz et al.           0.6011 [0.3914; 0.8108] < 0.0001  0.0878
## Omitting Kang et al.            0.5421 [0.3489; 0.7354] < 0.0001  0.0655
## Omitting Kuhlmann et al.        0.6055 [0.4048; 0.8063] < 0.0001  0.0746
## Omitting Lever Taylor et al.    0.5926 [0.3807; 0.8045] < 0.0001  0.0922
## Omitting Phang et al.           0.5840 [0.3707; 0.7972] < 0.0001  0.0933
## Omitting Rasanen et al.         0.5896 [0.3776; 0.8016] < 0.0001  0.0922
## Omitting Ratanasiripong         0.5829 [0.3719; 0.7939] < 0.0001  0.0900
## Omitting Shapiro et al.         0.5211 [0.3436; 0.6986] < 0.0001  0.0468
## Omitting Song & Lindquist       0.5798 [0.3661; 0.7934] < 0.0001  0.0934
## Omitting Warnecke et al.        0.5804 [0.3672; 0.7935] < 0.0001  0.0927
##                                                                         
## Pooled estimate                 0.5771 [0.3782; 0.7760] < 0.0001  0.0820
##                                     tau    I^2
## Omitting Call et al.             0.3003  64.2%
## Omitting Cavanagh et al.         0.3027  64.5%
## Omitting DanitzOrsillo           0.1882  48.1%
## Omitting de Vibe et al.          0.2716  57.6%
## Omitting Frazier et al.          0.3081  64.7%
## Omitting Frogeli et al.          0.3068  64.4%
## Omitting Gallego et al.          0.3007  63.8%
## Omitting Hazlett-Stevens & Oren  0.3076  64.8%
## Omitting Hintz et al.            0.2962  63.6%
## Omitting Kang et al.             0.2559  59.8%
## Omitting Kuhlmann et al.         0.2731  61.4%
## Omitting Lever Taylor et al.     0.3036  64.7%
## Omitting Phang et al.            0.3055  64.8%
## Omitting Rasanen et al.          0.3037  64.8%
## Omitting Ratanasiripong          0.3000  64.8%
## Omitting Shapiro et al.          0.2163  54.6%
## Omitting Song & Lindquist        0.3055  64.6%
## Omitting Warnecke et al.         0.3045  64.7%
##                                               
## Pooled estimate                  0.2863  62.6%
## 
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Hartung-Knapp adjustment for random effects model (df = 16, 17)

Le package dmetar permet d’obtenir une visualisation graphique sous la forme de forest plot de ces résultats, en les ayant au préalable ordonné. Il faut en premier calculer les valeurs influentes avec la fonction InfluentialAnalysis et demander ensuite de faire une representation graphique de la sortie de résultats avec la fonction plot (voir la Figure 11 et 12 pour les représentations des tailles d’effet et du \(I^2\)).

m.gen.inf<- InfluenceAnalysis(m.gen, random = TRUE)
## [===========================================================================] DONE
plot(m.gen.inf, "es")

plot(m.gen.inf, "i2")

Ainsi, sur la Figure qui montre l’évolution de la taille d’effet, l’étude qui est le plus haut est l’étude, qui après avoir été supprimée va entraîné la plus forte diminution de la taille d’effet, et celle qui est en bas, celle qui entraîne sa plus forte augmentation. Pour la Figure sur le \(I^2\), l’étude le plus haut indique celle qui va entraîner la plus forte diminution de l’hétérogénéité. Dans ces graphiques, la ligne verticale hachurée représente la taille d’effet originale et la zone verte représente l’intervalle de confiance sur la taille d’effet originale. Ainsi, on peut aisément identifier si une étude est clairement située en dehors de cet intervalle de confiance. Sans surprise, ce sont à nouveau les deux mêmes études qui ont le plus d’influence, et sont placées en haut du graphique.

La valeur chapeau

La valeur chapeau est la diagonale de la matrice d’influence. Il s’agit de l’écart entre les valeurs observées et les valeurs prédites. Cette mesure décrit l’influence que chaque mesure a sur les valeurs ajustées. Ainsi, la diagonale représente le levier qu’exercice chaque étude sur la taille d’effet moyenne.

On considère qu’une observation exerce une influence importante dès lors que sa valeur est supérieure à :

\[3\times \frac{p}{k}\text{ (5)}\]

Ainsi, une valeur avec être considérée comme influente si la valeur chapeau est supérieure à 0.1666667. On peut faire la représentation plus précise de l’évolution de cette mesure à l’aide de la fonction plot.

plot(inf.out, plotinf=7)

Figure 13. Graphique de la valeur chapeau.

Il faut noter ici que l’étude 3 reste identifiée comme influente, même si, sur la base du critère établi, elle ne devrait pas être considérée comme telle.

Le poids de l’étude

Nous avons déjà évoqué le poids des études dans le chapitre consacré aux effet fixes et aléatoires.

Nous nous limiterons donc en faire la présentation graphique.

plot(inf.out, plotinf=8)

Figure 14. Poids des études.

Cette mesure ne nécessite pas de description supplémentaire.

Le DFBETAS

Il existe une \(9^{ème}\) mesure qui n’est pas présentée par défaut, à savoir le DFBETAS. Cette mesure a plus de sens dans le contexte des méta-régresssions mais nous l’abordons dès à présent pour ne pas revenir sur la question des mesures d’influence plus loin.

Cette valeur indique les changements, en termes d’écart-types, qu’on observe sur le ou les coefficients estimés quand on compare le modèle après avoir ou non exclu une observation k.

Pour avoir cette mesure, il faut la demander explicitement avec la fonction dfbetas. On peut aussi la demander graphiquement avec la fonction plot, voir Figure 15 :

dfbetas(rma.out)
## 
##                        intrcpt 
## Call et al.             0.0365 
## Cavanagh et al.        -0.1955 
## DanitzOrsillo           0.9443 
## de Vibe et al.         -0.3321 
## Frazier et al.         -0.1795 
## Frogeli et al.         -0.0165 
## Gallego et al.          0.0547 
## Hazlett-Stevens & Oren -0.0855 
## Hintz et al.           -0.2508 
## Kang et al.             0.3910 
## Kuhlmann et al.        -0.3102 
## Lever Taylor et al.    -0.1597 
## Phang et al.           -0.0699 
## Rasanen et al.         -0.1282 
## Ratanasiripong         -0.0594 
## Shapiro et al.          0.6822 
## Song & Lindquist       -0.0272 
## Warnecke et al.        -0.0333
plot(inf.out, plotinf=F, plotdfbs =T)

Figure 15. Graphique des dfbetas.

Ainsi, il y a un dfbeta pour chaque predicteur de la méta-régression et les valeurs les plus larges indiquent la présence d’une valeur influente. Les valeurs supérieures à 1 constitue le critère utilisé pour identifier les valeurs comme influentes (Viechtbauer & Cheung, 2010).

Les mesures d’influence synthèse

Dès lors qu’une observation est considérée comme extrême et influente, elle peut affecter de manière négative les robustesse de la méta-anaalyse. Les règles pour savoir si une observation est considérée comme influente sont des heuristiques, mais comme toujours, il est nécessaire de comprendre comment une ou plusieurs études ont une influence sur vos réultats dans le contexte spécifique de votre étude. Ainsi, au-delà des heuristiques, cette analye réflexive est sans doute l’approche la plus adaptée pour savoir s’il faut ou non supprimer une étude du calcul de l’effet moyen pondéré (Harrer et al., 2021).

Ainsi, si on en reste à l’application par défaut des heuristiques définies, seule l’étude de DanitzOrsillo est considérée comme influente. Pourtant, dans tous les graphiques, on observe aussi que l’étude de Shapiro et al. forme un pic par rapport aux autres études, c’est-à-dire les études que nous avions identifiées avec le graphique de Baujat.

Ainsi, pour Harrer (2021), il y a un ensemble de preuves qui semblent indiquer qu’une grande part de l’hétérogénéité qui était initialement présente entre les études pourrait être causée par ces deux études.

Cependant, le problème avec cette logique est que nous ne pouvons pas fixer les critères pour identifier les valeurs influentes a priori si nous nous laissons guider par les données et on augmente dès lors le risque de se laisser guider par notre agenda, à savoir supprimer des études parce que leurs résultats ne nous plaisent pas.

L’analyse du graphique de Gosh

Une autre manière d’explorer l’hétérogénéité des données est d’utiliser une représentation graphique de l’hétérogénéité, appelé un graphique de GOSH (Olkin2012?).Pour ce type de graphique, on ajuste le même modèle de méta-analyse pour l’intégralité des sous-ensembles d’études que nous avons incluses. Ainsi, on va faire une méta-analyse avec toutes les études (k), une méta-analyse avec toutes les combinaisons de k-1 études, une méta-analyse avec toutes les combinaisons de k-2 études … (Viechtbauer & Cheung, 2010). Cela revient donc à estimer \(2^k-1\) modèles.

Cela signifie que, si le nombre d’études est important, le temps de calcul peut être considérablement long (Harrer et al., 2021). Ainsi, s’il y a 20 études, cela revient à estimer 1048575 modèles. C’est pourquoi le nombre de modèles estimés est limité à \(10^6\) modèles. Dans ce cas, c’est une sélection aléatoires des modèles possibles qui sera utilisée.

Si les tailles d’effet sont homogènes, le graphique va être approximativement symétrique et la distribution sera unimodale. Si la distribution est multimodale, cela suggère la présence d’une hétérogénéité, possiblement due à la présence de valeurs influentes ou de sous-groupes distincts d’études.

Une fois le modèle estimé, on peut montrer la dispersion des tailles d’effet sur l’axe des x et l’hétérogénéité entre les études sur l’axe des y. Ce qui est intéressant avec cette procédure est le fait que cela amène à pouvoir identifier des clusters d’études, dont chaque cluster a une taille d’effet différente. Il montre également l’étendue de l’hétérogénéité.

Un graphique de GOSH avec différents clusters indique qu’il y a plusieurs tailles d’effet de population, suggérant la pertinence de faire des analyses par sous-groupes.

Pour obtenir le graphique de GOSH, on utilise la fonction gosh sur la sortie de résultats de la méta-analyse (voir Figure 16. L’argument alpha indique la transparence des points. Comme il y a beaucoup de points, on n’est pas en mesure de détecter les points qui s’empilent les uns sur les autres. Le fait d’utiliser une transparence permet de contourner ce souci.

res.gosh<-gosh(rma.out)
res.gosh
## 
## Model fits attempted: 262143
## Model fits succeeded: 261709
## 
##             mean    min      q1  median      q3     max 
## k         9.0000 1.0000  8.0000  9.0000 10.0000 18.0000 
## QE       21.4759 0.0000 13.9468 21.5322 28.5681 45.5026 
## I2       56.1321 0.0000 42.6133 64.8774 75.2550 95.6523 
## H2        3.0552 1.0000  1.7426  2.8472  4.0412 23.0007 
## tau2      0.0979 0.0000  0.0336  0.0869  0.1481  1.3453 
##                                                         
## tau       0.2794 0.0000  0.1832  0.2948  0.3848  1.1599 
## estimate  0.5853 0.1036  0.5077  0.5821  0.6580  1.7912
plot(res.gosh, alpha = 0.01)

Figure 16. Graphique de Gosh.

On observe que la distribution du \(I^2\) est asymétrique et bimodale. Il semble donc qu’il existe des combinaisons d’études pour lesquelles l’hétérogénéité est bien plus faible. De manière interessant, la forme d’étoile filante que prend le graphique indique aussi que les tailles d’effet sont bien plus faibles lorsque l’homogénéité est plus faible.

Pour identifier les études qui sont à l’origine de cette asymétrie, on peut utiliser la fonction gosh.dianostics issu du package dmetar (Harrer et al., 2019b). Cette fonction utilise un algorithme d’apprentissage machine non supervisé pour détecter les clusters dans les données du graphique de GOSH. Il va détecter automatiquement les études qui contribuent le plus à chacun des clusters. Si une étude est systématiquement, ou presque systématiquement présente dans un cluster alors que ce n’est pas le cas des autres, cela indique que cette étude (ou cette combinaison d’études) pourrait être la cause de l’hétérogénéité.

Les packages requis pour la fonction gosh.dianostics.

La fonction gosh.dianostics requiert d’autres packages qu’il faut avoir installé si ce n’est pas encore le cas et chargé. Ces packages sont ‘ggplot2’ (ggplot2?), ‘gridExtra’ (gridExtra?), ‘mclust’ (mclust?) et ‘fpc’ (fpc?).

Si ces packages ne sont pas installés, il faut utiliser les lignes de commande suivantes :

install.packages("ggplot2")
install.packages("gridExtra")
install.packages("mclust")
install.packages("fpc")

Pour charger ces packages, il faut utiliser les lignes de commande suivantes :

require("ggplot2")
require("gridExtra")
require("mclust")
require("fpc")

La fonction goh.diagnostic>2 va utiliser différents algorithmes pour détecter des clusters : la méthode basée sur le k-means (Hartigan1979?), la méthode DBSCAN (Schubert2017?) et les modèles gaussiens mixtes (Fraley2002?).

Il est possible de spécifier un certain nombre de paramètres pour chacun de ces algorithmes. Harrer et . (2019b) suggère, pour la méthode k-means, de rechercher 2 clusters dans cet exemple. Il faut cependant noter que les algorithmes ne sont pas optimaux dans tous les situations et les résultats de cette fonction doivent être pris avec précaution si vous n’êtes pas en mesure de customiser les paramètres.

res.gosh.diag <- gosh.diagnostics(res.gosh, 
                                  km.params = list(centers = 2), # on indique le nombre de centres sur la méthode k-means
                                  db.params = list(eps = 0.08, 
                                                   MinPts = 50))
res.gosh.diag

La sortie de résultats indique le nombre de clusters que chaque méthode a détecté et l’identification de valeurs potentiellement aberrantes/influentes. C’est tout à fait normal que le nombre de clusters ne soit pas le même pour les différentes approches. Ce qui est intéressant, c’est d’identifier que, parmi les valeurs influentes, il apparaît que les études 3 et 16, les mêmes que nous avions détectées précédemment sont également détectées ici comme outliers.

plot<-plot+1

Pour avoir une représentation graphique de ces résultats, on peut utiliser la fonction plot (voir Figure 17).

plot(res.gosh.diag)

Figure 17. Résultats des méthodes de clustering sur les résultats de la figure de GOSH

Les trois premiers graphiques présente les solution de clustering identifiées par chaque algorithme et le niveau de déséquilibre associée à chacune des étude. Sur la base de ces informations, une distance de Cook est calculée afin de déterminer si une étude est suscepptible d’avoir un grand impact sur la détection du cluster, ce qui suggérerait qu’elle est influente.

Les autres graphiques sont à nouveau les graphiques de GOSH.

Une étude qui mérite notre attention est celle de Vibe qui pourrait être considérée comme influente aussi. Ceci s’explique parce que la taille d’échantillon est particulièrement grande, entraînant une mesure très précise de la taille d’effet et donc un poids plus conséquent sur les résultats.

Que fait quand on a détecté des valeurs influentes

Quand des valeurs influentes ont été détectées, il est raisonnable de les supprimer. Il faut donc refaire la méta-analyse en excluant ces études. Néanmoins, par souci de transparence, le mieux quand on présente les résultats est d’indiquer les résultats avec et sans ces valeurs influentes.Pour cela, on va utiliser l’argument subset de la fonction rma.

rma.out.sans.infl<-rma(yi = TE, sei = seTE, data =ThirdWave, method="REML" ,test="knha", slab=Author,
             subset = c(1:2,4:15,17:18))

rma.out.sans.infl
## 
## Random-Effects Model (k = 16; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0139 (SE = 0.0200)
## tau (square root of estimated tau^2 value):      0.1180
## I^2 (total heterogeneity / total variability):   25.00%
## H^2 (total variability / sampling variability):  1.33
## 
## Test for Heterogeneity:
## Q(df = 15) = 19.9498, p-val = 0.1739
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub      
##   0.4528  0.0597  7.5897  15  <.0001  0.3257  0.5800  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En supprimant les deux études identifiées comme outlier, l’\(I^2\) tombe à 25% (et même à 0 si on supprimait l’étude de Vibe). Que ce soit dans un cas ou dans l’autre, la taille d’effet est sensiblement la même : aux environ de 0.45. Cette valeur est plus petite que la taille d’effet initialement identifiée, même si on reste sur le même ordre de grandeur. Donc, même avec les valeurs influentes, nos résultats n’étaient pas démesurément biaisées par ces valeurs.

Lorsqu’on reporte les résultats, il peut être utile pour le lecteur d’être en mesure de comprendre quels changements ont eu la suppression des valeurs influentes, c’est pourquoi Harrer et al. (Harrer et al., 2021) recommandent de faire un tableau avec les changements sur les mesures d’hétérogénéité, tout en spécifiant les études qui ont été supprimées des analyses.

Les résultats peuvent dès lors être présentés de la manière suivante :

Nous avons mené une méta-analyse sur les TCC de \(3^{ème}\) génération sur le stress perçu chez les étudiants par rapport à un groupe de contrôle inactif. L’analyse de l’hétérogénéité suggèrait la présence de valeurs influentes ou de plusieurs clusters. Une analyse plus détaillées des valeurs influentes indiquait que 2 études avaient une influence considérable sur les résultats. Les résultats sont résumés dans le Tableau 4.

require(flextable)
require(officer)
IC95.1 <-paste0( round(rma.out.sans.infl$ci.lb,3), " − ",round(rma.out.sans.infl$ci.ub,3) )
        
IC95.2<-paste0( round(rma.out.sans.infl$ci.lb,3), " − ",round(rma.out.sans.infl$ci.ub,3) )
pval1<-ifelse(round(rma.out$pval, 3)==0, "<.001", rma.out$pval)
pval2<-ifelse(round(rma.out.sans.infl$pval, 3)==0, "<.001", rma.out$pval)

PI1<-paste0( round(predict(rma.out)$pi.lb,2), " − ",
             round(predict(rma.out)$pi.ub ,2))
pI2<-paste0(round(predict(rma.out.sans.infl)$pi.lb,2), 
            " − ",
            round(predict(rma.out.sans.infl)$pi.ub ,2))


I.CI1<-paste0(round(confint(rma.out)$random[3,2],0), "% − ", 
       round(confint(rma.out)$random[3,3],0), "%" )
I.CI2<-paste0(round(confint(rma.out.sans.infl)$random[3,2],0), "% − ", 
       round(confint(rma.out.sans.infl)$random[3,3],0), "%" )

out<- rma.out$slab[which(rma.out.sans.infl$subset==F)] 
out<-paste0(out, collapse= " and ")
footnote<-paste0("Studies from ", out, " have been considered as influential.")
resultats<-data.frame(Analyse = c("Analyse principale", "Analyse sans les valeurs influentes"),
                      g = c(round(rma.out$b,3), round(rma.out.sans.infl$b,3)),
                      IC95 =c(IC95.1,IC95.2),
                      p = c(pval1, pval2),
                      PI95= c(PI1,pI2),
                        I2 = c(round(rma.out$I2,2), round(rma.out.sans.infl$I2,2)),
                      I.CI95 = c(I.CI1, I.CI2)   )
myft<-flextable( resultats)
myft<-set_table_properties(myft, layout = "autofit")
myft<-set_header_labels(x=myft, values = c(" ", "g", "CI 95%",
                                           "p", "PI 95%",
                                           "I²", "I² CI 95%"))
myft <- add_footer_lines(myft, footnote)
myft<-  style(myft,
              pr_p = fp_par(
                text.align = "center", padding = 1))
myft<- set_caption(myft, caption = "Résultats de la méta-analyse sur les TCC de 3ème génération")                
myft
Résultats de la méta-analyse sur les TCC de 3ème génération

g

CI 95%

p

PI 95%

I²

I² CI 95%

Analyse principale

0.577

0.326 − 0.58

<.001

-0.06 − 1.21

64.63

40% − 89%

Analyse sans les valeurs influentes

0.453

0.326 − 0.58

<.001

0.17 − 0.73

25.00

0% − 71%

Studies from DanitzOrsillo and Shapiro et al. have been considered as influential.

Astuce

Le tableau ci-dessus est directement utilisable pour les résultats. Il n’y a que deux contraintes à son utilisation par copier coller :

  1. Les sorties de résultats de la fonction rmarma.out et rma.out.sans.infl
  2. Les packages ‘flextable’ (Gohel & Skintzos, 2024) et officer (Gohel & Moog, 2024) doivent être installés et chargés.

Il faut néanmoins avoir conscience que cela reste de votre responsabilité de vérifier que les valeurs reportées sont correctes.

Exercice 1

Reprenez l’exercice sur l’efficacité du vaccin BCG présenté au chapitre précédent (pour rappel, les données sont directement disponibles dans le package “metafor”). Calculez différents indices permettant d’identifier la présence de valeurs influentes et/ou de clusters. Présentez les résultats globaux, et si cela est pertinent sans les valeurs influentes (Pour cet exercice, vous pouvez ignorer le forest plot).

Cliquez pour la solution

Les données sont présentées ci-dessous (voir Tableau 5.

Données du vaccin BCG

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

On commence par calculer les tailles d’effets.

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

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

rma.out <- rma(yi, vi, data=dat, method="PM", digits=2, slab = author)
rma.out
## 
## Random-Effects Model (k = 13; tau^2 estimator: PM)
## 
## 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

On observe ici que le \(I^2\) vaut 92.15%, suggérant la pertinence d’aller explorer plus en détail s’il y a des valeurs influentes.

inf <- influence(rma.out)
inf
## 
##                      rstudent dffits cook.d cov.r tau2.del QE.del  hat weight 
## Aronson                 -0.23  -0.04   0.00  1.15     0.38 162.45 0.05   5.00 
## Ferguson & Simes        -1.32  -0.35   0.12  1.00     0.31 155.83 0.06   6.35 
## Rosenthal et al.1       -0.74  -0.16   0.02  1.09     0.36 161.07 0.05   4.50 
## Hart & Sutherland       -1.29  -0.44   0.18  1.04     0.32 106.98 0.10   9.65 
## Frimodt-Moller et al     0.87   0.28   0.08  1.12     0.35 162.23 0.09   8.87 
## Stein & Aronson         -0.36  -0.10   0.01  1.20     0.38 129.61 0.10   9.94 
## Vandiviere et al        -1.24  -0.32   0.10  1.01     0.32 156.80 0.06   6.14 
## TPT Madras               1.43   0.47   0.20  1.02     0.30  72.72 0.10  10.11 
## Coetzee & Berjak         0.44   0.15   0.02  1.18     0.37 163.14 0.09   8.76 
## Rosenthal et al.2       -1.07  -0.33   0.10  1.08     0.34 150.52 0.08   8.37 
## Comstock et al.1         0.70   0.24   0.06  1.16     0.36 162.32 0.10   9.86 
## Comstock & Webster       1.32   0.25   0.06  0.97     0.31 161.70 0.04   3.99 
## Comstock et al.2         1.21   0.36   0.12  1.05     0.32 160.65 0.08   8.45 
##                       dfbs inf 
## Aronson              -0.04     
## Ferguson & Simes     -0.35     
## Rosenthal et al.1    -0.15     
## Hart & Sutherland    -0.44     
## Frimodt-Moller et al  0.28     
## Stein & Aronson      -0.10     
## Vandiviere et al     -0.32     
## TPT Madras            0.47     
## Coetzee & Berjak      0.15     
## Rosenthal et al.2    -0.33     
## Comstock et al.1      0.24     
## Comstock & Webster    0.25     
## Comstock et al.2      0.36
### plot the values
plot(inf)

Aucune étude ne semble être influente en l’occurrence. Nous utilisons un seuil de 0.83 pour le DFFIT, 0.45 pour la distance de Cook et 0.23 pour la valeur chapeau. On remarque néanmoins que, si on supprime l’étude 8, cela diminue considérablement l’hétérogénéité.

Nous pouvons continuer avec un graphique de Baujat et celui de GOSH.

baujat(rma.out, symbol="slab")

Le graphique de Baujat ne permet pas de clairement identifier des valeurs influentes, bien que l’étude de Hart et Sutherlan et celle de TPT Madras semblent s’éloigner un peu plus. L’étude de TPT Madras avait déjà été identifiée comme étant celle qui avait le plus gros impact sur l’hétérogénéité.

gosh.out<-gosh(rma.out)

gosh.out
## 
## Model fits attempted: 8191
## Model fits succeeded: 8191
## 
##           mean   min    q1 median    q3    max 
## k         6.50  1.00  5.00   7.00  8.00  13.00 
## QE       62.18  0.00 24.39  47.43 98.42 163.16 
## I2       85.25  0.00 83.43  89.55 93.06  98.87 
## H2       10.71  1.00  6.03   9.57 14.42  88.66 
## tau2      0.34  0.00  0.27   0.34  0.42   1.86 
##                                                
## tau       0.57  0.00  0.52   0.59  0.65   1.36 
## estimate -0.75 -1.67 -0.87  -0.75 -0.62   0.45

Etant donné que nous avions identifié l’étude 8 comme potentielle valeur influente, nous pouvons regarder la distribution de la taille d’effet de l’hétérogénéité lorsque cette étude est présente (en bleu) ou non (en rouge).

plot(gosh.out, out=8, breaks=100)

On se rend compte qu’avec ou sans cette étude, la taille d’effet estimé est sensiblement la même, la distribution est symétrique et unimodale. En revanche, supprimer cette étude diminue conidérablement l’hétérogénéité des résultats.

Pour l’exemple pédagogique, nous allons donc continuer en réajustant la méta-analyse sans cette valeur et nous présenterons ensuite les résultats.

rma.out.sans.infl<- rma(yi, vi, data=dat, method="PM", digits=2, slab = author,
             subset = c(1:7, 9:13))
table<-table+1

Nous avons mené une méta-analyse sur l’efficacité du vaccin BCG contre la tuberculose. L’analyse de l’hétérogénéité suggèrait la présence de valeurs influentes ou de plusieurs clusters. Une analyse plus détaillées des valeurs influentes indiquait que l’étude augmentait considérablement l’hétérogénéité. Les résultats sont résumés dans le Tableau 6.

require(flextable)
require(officer)
IC95.1 <-paste0( round(rma.out.sans.infl$ci.lb,3), " − ",round(rma.out.sans.infl$ci.ub,3) )
        
IC95.2<-paste0( round(rma.out.sans.infl$ci.lb,3), " − ",round(rma.out.sans.infl$ci.ub,3) )
pval1<-ifelse(round(rma.out$pval, 3)==0, "<.001", rma.out$pval)
pval2<-ifelse(round(rma.out.sans.infl$pval, 3)==0, "<.001", rma.out$pval)

PI1<-paste0( round(predict(rma.out)$pi.lb,2), " − ",
             round(predict(rma.out)$pi.ub ,2))
pI2<-paste0(round(predict(rma.out.sans.infl)$pi.lb,2), 
            " − ",
            round(predict(rma.out.sans.infl)$pi.ub ,2))


I.CI1<-paste0(round(confint(rma.out)$random[3,2],0), "% − ", 
       round(confint(rma.out)$random[3,3],0), "%" )
I.CI2<-paste0(round(confint(rma.out.sans.infl)$random[3,2],0), "% − ", 
       round(confint(rma.out.sans.infl)$random[3,3],0), "%" )

out<- rma.out$slab[which(rma.out.sans.infl$subset==F)] 
out<-paste0(out, collapse= " and ")
footnote<-paste0("Studies from ", out, " have been considered as influential.")
resultats<-data.frame(Analyse = c("Analyse principale", "Analyse sans les valeurs influentes"),
                      g = c(round(rma.out$b,3), round(rma.out.sans.infl$b,3)),
                      IC95 =c(IC95.1,IC95.2),
                      p = c(pval1, pval2),
                      PI95= c(PI1,pI2),
                        I2 = c(round(rma.out$I2,2), round(rma.out.sans.infl$I2,2)),
                      I.CI95 = c(I.CI1, I.CI2)   )
myft<-flextable( resultats)
myft<-set_table_properties(myft, layout = "autofit")
myft<-set_header_labels(x=myft, values = c(" ", "g", "CI 95%",
                                           "p", "PI 95%",
                                           "I²", "I² CI 95%"))
myft <- add_footer_lines(myft, footnote)
myft<-  style(myft,
              pr_p = fp_par(
                text.align = "center", padding = 1))
myft<- set_caption(myft, caption = "Résultats de la méta-analyse sur les TCC de 3ème génération")                
myft
Résultats de la méta-analyse sur les TCC de 3ème génération

g

CI 95%

p

PI 95%

I²

I² CI 95%

Analyse principale

-0.745

-1.198 − -0.46

<.001

-1.95 − 0.46

92.15

82% − 98%

Analyse sans les valeurs influentes

-0.829

-1.198 − -0.46

<.001

-1.97 − 0.31

87.16

69% − 96%

Studies from TPT Madras have been considered as influential.

Références

Balduzzi, S., Rücker, G., & Schwarzer, G. (2019). How to perform a meta-analysis with R: A practical tutorial. Evidence-Based Mental Health, 22, 153–160.
Baujat, B., Mahé, C., Pignon, J.-P., & Hill, C. (2002). A graphical method for exploring heterogeneity in meta-analyses: Application to a meta-analysis of 65 trials. Stat. Med., 21(18), 2641–2652.
Gohel, D., & Moog, S. (2024). Officer: Manipulation of microsoft word and PowerPoint documents. https://CRAN.R-project.org/package=officer
Gohel, D., & Skintzos, P. (2024). Flextable: Functions for tabular reporting. https://CRAN.R-project.org/package=flextable
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. (2019b). Dmetar: Companion r package for the guide ’doing meta-analysis in r’. http://dmetar.protectlab.org/
Harrer, M., Cuijpers, P., Furukawa, T., & Ebert, D. D. (2019a). Dmetar: Companion r package for the guide ’doing meta-analysis in r’. http://dmetar.protectlab.org/
Knapp, G., & Hartung, J. (2003). Improved tests for a random effects meta‐regression with a single covariate. Statistics in Medicine, 22(17), 2693–2710. https://doi.org/10.1002/sim.1482
Viechtbauer, W. (2020). Model checking in meta-analysis. In C. H. Schmid, T. Stijnen, & I. White (Eds.), Handbook of meta-analysis (pp. 219–254). Chapman & Hall/CRC.
Viechtbauer, W., & Cheung, M. W.-L. (2010). Outlier and influence diagnostics for meta-analysis. Research Synthesis Methods, 1(2), 112–125. https://doi.org/10.1002/jrsm.11

  1. il est à noter que les calculs ont été légèrement simplifiés ici pour considérer qu’il s’agissait d’un modèle à effets fixes↩︎

  2. la fonction semble avoir un bug au moment de la rédaction de ces lignes mais l’auteur du package a été contacté et on peut espérer que la fonction soit à nouveau fonctionnelle prochainement.↩︎