TD Meta-Analyse

\(~\)

1 INTRODUCTION

1.1 Rappel d’une démarche de Méta-analyse

1. Définir clairement la question de recherche, systèmes étudiées, zone considérée
2. Recherche exhaustive de la litérature
3. Création de la base de données, et analyse de la qualité des études
4. Analyse quantitative des résultats

-> Ce TD va se focaliser sur la partie analyse quantitative

1.2 Analyse quantitative des résultats

Selon Philibert et al., 2012 (voir aussi Beillouin et al., 2019), l’analyse quantitative se décompose en plusieurs parties:

1. Analyse des résultats, et de la variabilité inter et intra etudes en prenant en compte (weighting) la qualiét et/ou précision des études
2. analyse de sensibilité des résultats aux modification du jeu de données
3. Analyse du biais de publication

2 DEBUT DU TD

2.1 On charge la liste des packages

#install.packages("KenSyn")
library(KenSyn)   # <- pour le jeu de données
library(ggplot2)  # <- pour les représentations graphiques
library(psych)    # <- pour les représentations graphiques
library(gmodels)  # Compléments
library(metafor)  # Package pour la méta-analyse
library(tidyverse) # Pour coder plus facilement

2.2 On va travailler à partir d’un example

Les données ont été récupérées à partir de la littérature pour une méta-analyse visant à comparer le rendement des cultures en système biologique au système conventionnel effectuée par Seufert et al (2012). Il contient 65 études. Chaque étude contient des valeurs de rendement des systèmes biologique et conventionnel comparées dans le cadre d’expériences sur le terrain. Les études couvrent différents lieux géographiques, années et cultures.

Data <- organic

2.3 Analyse graphique et description des données

2.3.1 Description à plat des données

library(pastecs)
describe(Data)

2.3.2 Rendement bio vs. conventionnels

p<-ggplot(Data, aes(Yield_conv, Yield_org,, color=Crop_species)) +
  geom_point() +
  theme_classic()
ggExtra::ggMarginal(p, type = "histogram")

2.3.3 Histogramme des ratios de rendements

ggplot(Data, aes(lnR)) +
  geom_vline(aes(xintercept=0), linetype=2, color= "gray50")+
  geom_histogram(color="black", fill="gray80")+
  geom_density(alpha=.2, fill="#FF6666") 

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

2.3.4 Forest plot des ratios de rendements par étude

p<-ggplot(Data) +
  geom_hline(aes(yintercept=0), linetype=2, color= "gray50")+
  geom_point(aes(reorder(Study, lnR),lnR, size= N_conv,color=Crop_species), fill="gray80")+
  geom_errorbar(aes(reorder(Study, lnR),ymin = lnR-1.96*Var_lnR, ymax = lnR+Var_lnR), width = 0.2)+
  coord_flip()+
  geom_text(aes(Study,-1.5, label= paste("n = ", N_conv)))
p

2.3.5 Quelles co-variables semblent influencer le ratio de rendement?

Data_NUM <- Data %>% dplyr::select_if(is.numeric)
pairs.panels(Data_NUM, 
             method = "pearson", # correlation method
             hist.col = "#00AFBB",
             density = TRUE,  #  show density plots
             ellipses = FALSE # NOT show correlation ellipses
             )

library(ranger) Data_RF <- Data %>% dplyr:: select(-Study,-Yield_conv,-SD_conv,-N_conv,-Yield_org,-SD_org,-N_org, -Yield_unit, -yi, -vi,-Var_lnR) data=organic library(ranger) mod<-ranger(lnR~., data=Data_RF,importance = “permutation”) library(vip) vip(mod)

2.4 Calculs des effect-sizes

2.4.1 Fixed effect models

On va utiliser la fonction metafor. Une description complete de la fonction est accéssible en tapant dans la console ?metafor ou vignette(“metafor”)

ma_model_1 <- rma(lnR, Var_lnR, method="FE", data = organic)
summary(ma_model_1)

## 
## Fixed-Effects Model (k = 65)
## 
##    logLik   deviance        AIC        BIC       AICc 
## -922.5790  1982.3103  1847.1579  1849.3323  1847.2214   
## 
## Test for Heterogeneity:
## Q(df = 64) = 1982.3103, p-val < .0001
## 
## Model Results:
## 
## estimate      se      zval    pval    ci.lb    ci.ub 
##  -0.2528  0.0082  -30.7315  <.0001  -0.2690  -0.2367  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

-> Interpretation des sorties du modèle

1. AIC, BIC, …. : Indicateur de la qualité d’ajustement du modèle. Peux servir à comparer d’autres modèles, mais sur la même base de données

2. Test for Heterogeneity:
   Q = ∑ (Wi(Yi-θ)^2, with Wi : Weight, Yi: observed effect-size, θ: estimated effect size Q is the weighted deviations about the summary effect size. Larger values of Q reflect greater between-study heterogeneity. Q-test has low power (<0.80) when the number of studies and/or sample sizes is small.

3. p-value:
   pchisq(Q, df = df, lower = FALSE)
   how likely under null hypothesis (H0:θ1= θ2= …)
   p-val < .0001 -> rejet H0.
   

4. Model results estimate: -0.2528. le log ratio de rendement BIO/CONV est -0.24 <-> le ratio de rendement est exp(-0.25)= 0.7766 <-> le Bio produit 100-77= 22.4% de moins que le conventionnel A quel point ce résultat est précis? on regarde les confidence intervals la moyenne de perte de Rdt en Bio se situe entre 21.02- et 23.58% par rapport au conventionnel

-> Représentation graphique

directement avec la fonction metafor

forest(ma_model_1)

ou avec le package ggplot2

RES<-summary(ma_model_1)
Data["MA_fixed", "lnR"]<- RES$beta
Data["MA_fixed", "Var_lnR"]<- RES$se^2
Data["MA_fixed", "Study"]<- "Fixed_effect summary"

p<-ggplot(Data) +
  geom_hline(aes(yintercept=0), linetype=2, color= "gray50")+
  geom_point(aes(reorder(Study, lnR),lnR, size= N_conv,color=Crop_species), fill="gray80")+
  geom_errorbar(aes(reorder(Study, lnR),ymin = lnR-1.96*Var_lnR, ymax = lnR+Var_lnR), width = 0.2)+
  coord_flip()+
  geom_text(aes(Study,-1.5, label= paste("n = ", N_conv)))
p

## Warning: Removed 1 rows containing missing values (geom_point).

2.4.2 Random effect models

Il existe différentes façon d’écrire les modèles mixtes (avec metafor)

  ma_model_3<-rma.mv(lnR, Var_lnR, random = ~ 1 | Study, data=organic)
  summary(ma_model_3)

## 
## Multivariate Meta-Analysis Model (k = 65; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -20.0769   40.1539   44.1539   48.4717   44.3506   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed  factor 
## sigma^2    0.0724  0.2690     65     no   Study 
## 
## Test for Heterogeneity:
## Q(df = 64) = 1982.3103, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub 
##  -0.2429  0.0398  -6.1004  <.0001  -0.3209  -0.1649  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  ma_model_2<- rma(lnR, Var_lnR, data=organic)
  summary(ma_model_2)

## 
## Random-Effects Model (k = 65; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -20.0769   40.1539   44.1539   48.4717   44.3506   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0724 (SE = 0.0171)
## tau (square root of estimated tau^2 value):      0.2690
## I^2 (total heterogeneity / total variability):   93.47%
## H^2 (total variability / sampling variability):  15.31
## 
## Test for Heterogeneity:
## Q(df = 64) = 1982.3103, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub 
##  -0.2429  0.0398  -6.1004  <.0001  -0.3209  -0.1649  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

-> Interpretation des sorties du modèle

1. Les résultats indiquent que le rapport logarithmique moyen estimé est égal à μˆ = -0,2429 (IC 95 % : -0,3203 à -0,1655). Pour faciliter l’interprétation, il peut être utile de retransformer ces valeurs à l’échelle des ratios par exponentiation (c.-à-d. exp(μˆ) = 0,78 avec IC à 95 % : 0,73 à 0,84). L’hypothèse nulle H0 : μ = 0 peut être clairement rejetée (z = -6,1503, p < 0,0001).

2. Intéressant de remarquer que ces ratios sont un peu différent de ceux de l’article Seufert et al., 2012. Les auteurs se basent sur 316 points, tandisque nous en avons 346. Ils ont du sélectioner les données et pondérer les études? (voir plus bas)

3. La statistique I2 estime (en pourcentage) quelle part de la variabilité totale des estimations de la valeur de l’effet (qui se compose de l’hétérogénéité et de la variabilité d’échantillonnage) peut être attribuée à l’hétérogénéité des effets réels (τˆ2 = 0 implique donc I2 = 0%).

4. La statistique H2 est le rapport entre la variabilité totale des résultats observés et la variabilité de l’échantillonnage (τˆ2 = 0 implique donc H2 = 1).

5. Il est important de comprendre, cependant, que τˆ2, I2 et H2 sont souvent estimés de façon imprécise, surtout lorsque le nombre d’études est faible. On peut calculer les intervalles des confiances sur ces indicateurs

  confint(ma_model_2)

## 
##        estimate   ci.lb   ci.ub 
## tau^2    0.0724  0.0437  0.1136 
## tau      0.2690  0.2091  0.3370 
## I^2(%)  93.4694 89.6327 95.7382 
## H^2     15.3125  9.6457 23.4641

2.4.3 Random effect models with weights

Pour ajuster les poids des études dans le modèle, on utilise la fonction weights de metafor

  ma_model_4<- rma(lnR, Var_lnR, data=organic, weights= 1/Var_lnR)
  summary(ma_model_4)

## 
## Random-Effects Model (k = 65; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -20.1081   40.2162   44.2162   48.5340   44.4129   
## 
## tau^2 (estimated amount of total heterogeneity): 0.0724 (SE = 0.0171)
## tau (square root of estimated tau^2 value):      0.2690
## I^2 (total heterogeneity / total variability):   93.47%
## H^2 (total variability / sampling variability):  15.31
## 
## Test for Heterogeneity:
## Q(df = 64) = 1982.3103, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub 
##  -0.2528  0.1021  -2.4764  0.0133  -0.4529  -0.0527  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

-> Interpretation des sorties du modèle

Avec le modèle 2, l’estimateur le plus efficace de mu est donné par
mu-hat = somme(wi * yi) / somme(wi)
où wi = 1/(tau^2 + vi), c’est-à-dire l’inverse des variances marginales

Dans certains cas, il se peut que nous voulions nous écarter de l’utilisation de ces poids. Supposons, par exemple, que l’on s’inquiète de la précision des études. . Si tau^2 est grand, alors mu-hat est presque comme une moyenne non pondérée et les études plus petites et plus grandes reçoivent presque le même poids dans le calcul mu-hat. Cela peut conduire à un biais substantiel dans le mu-hat. Au lieu de cela, nous pouvons estimer mu avec wi = 1/vi, de sorte que les petites études reçoivent peu de poids et les grandes études plus de poids. Par conséquent, le chapeau mu peut avoir beaucoup moins de biais.

2.4.4 Random effect models with moderaors

Pour cela on utilise l’argument mods de la fonction metafor
exemple en utilisant Crop_type, comme dans Seufert et al., 2012, Figure 1a.

ma_model_5 <- rma(lnR, Var_lnR, weights = 1/ Var_lnR, mods= ~Crop_type-1, data = organic)
summary(ma_model_5)

## 
## Mixed-Effects Model (k = 65; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc 
## -35.0465   70.0930   94.0930  117.9608  101.7028   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0809 (SE = 0.0204)
## tau (square root of estimated tau^2 value):             0.2845
## I^2 (residual heterogeneity / unaccounted variability): 93.26%
## H^2 (unaccounted variability / sampling variability):   14.83
## 
## Test for Residual Heterogeneity:
## QE(df = 54) = 1180.7318, p-val < .0001
## 
## Test of Moderators (coefficients 1:11):
## QM(df = 11) = 28.2389, p-val = 0.0030
## 
## Model Results:
## 
##                          estimate      se     zval    pval    ci.lb 
## Crop_typecereals          -0.1114  0.1622  -0.6867  0.4923  -0.4293 
## Crop_typefiber            -0.4032  0.2804  -1.4376  0.1506  -0.9528 
## Crop_typeforage            0.0436  0.2957   0.1474  0.8828  -0.5360 
## Crop_typefruits           -0.3208  0.1941  -1.6530  0.0983  -0.7012 
## Crop_typeoil crops        -0.1187  0.2164  -0.5484  0.5834  -0.5429 
## Crop_typeother            -0.3419  0.2333  -1.4655  0.1428  -0.7991 
## Crop_typepulses           -0.4518  0.5815  -0.7770  0.4372  -1.5915 
## Crop_typeroots & tubers   -0.4865  0.2916  -1.6683  0.0953  -1.0582 
## Crop_typesugar crops      -0.2930  0.4055  -0.7225  0.4700  -1.0878 
## Crop_typetree nuts        -0.0484  0.3560  -0.1360  0.8919  -0.7462 
## Crop_typevegetables       -0.6748  0.1658  -4.0707  <.0001  -0.9997 
##                            ci.ub 
## Crop_typecereals          0.2065      
## Crop_typefiber            0.1465      
## Crop_typeforage           0.6232      
## Crop_typefruits           0.0596    . 
## Crop_typeoil crops        0.3055      
## Crop_typeother            0.1154      
## Crop_typepulses           0.6879      
## Crop_typeroots & tubers   0.0851    . 
## Crop_typesugar crops      0.5018      
## Crop_typetree nuts        0.6494      
## Crop_typevegetables      -0.3499  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

-> Interpretation des sorties du modèle

1. Nous n’expliquons pas une large part de la variabilité (cf. comparaison τˆ2 dans les modèles 5 et 4). Probablement parce que le nombre d’individus dans certains groupes est faible -> mauvaise qualité des estimations.

2. Test des modérateurs :
   On peut rejeter H0 : β1 = β2 = 0 basé sur le test (QM = 28.2389, df = 11, p < 0.0030) : Les légumes de type crop_type semblent avoir une influence significative sur le rapport logarithmique (c.-à-d., pour H0 : β1 = 0, on trouve z=–4.0707 <0.01).

3. Test d’hétérogénéité résiduelle :
   Le test d’hétérogénéité résiduelle est significatif (QE = 1180, df = 54, p < 0,0001), ce qui indique peut-être que d’autres modérateurs non considérés dans le modèle influencent le rapport logarithmique.

-> représentation graphique

RES<-summary(ma_model_5)
SUMMARY    <- data.frame(b=exp(RES$beta))
SUMMARY$se<- RES$se
SUMMARY$low_CI<- exp(RES$b-RES$se)
SUMMARY$up_CI<- exp(RES$b+RES$se)
SUMMARY$Crop_type<- row.names(SUMMARY)

ggplot(data=SUMMARY)+
  geom_point(aes(reorder(Crop_type,b), b))+
  geom_errorbar(aes(Crop_type, ymin= low_CI, ymax=up_CI))+
  geom_hline(aes(yintercept=1), color= "gray80", linetype=2)+
  coord_flip()

2.5 Analyse de sensibilité des résultats des modèles

Il existe de nombreuses méthodes pour analyser la sensibilité des résultats aux études inclues dans la base de données. En pratiques ces méthodes ne sont que rarement mise en oeuvre dans les méta-analyses.

2.5.1 Diagnostic par influence

On va utiliser la fonction influence qui permet de réaliser rapidement une série de diagnostic

inf<-influence(ma_model_5)
plot(inf)

-> Interprétation des sorties

1. dffits: The DIFFITS value of a study indicates in standard deviations how much the predicted pooled effect changes after excluding this study.

2. cook’s distance: même idée, mais indicateur un peu différent.

3. cov.r: déterminant de la matrice de variance-covariance des estimations des paramètres lorsque l’étude est supprimée, divisé par le déterminant de la matrice de variance-covariance des estimations des paramètres lorsque l’ensemble des données est considéré. Il est important de noter que les valeurs de cov.r < 1 indiquent que la suppression de l’étude conduira à une estimation plus précise de la taille de l’effet (c.-à-d. moins hétérogénéité).

4. autre variables: voir Viechtbauer, W., & Cheung, M. W.-L. (2010). Outlier and influence diagnostics for meta-analysis. Research Synthesis Methods, 1, 112–125

5. L’élimination de ces études réduirait considérablement l’hétérogénéité et augmenterait la précision du résultat moyen estimé (c.-à-d. la corrélation transformée r-to-z.) du modèle à effets aléatoires. Cependant, au lieu de simplement supprimer ces études, on devrait les examiner en détail pour déterminer la raison de leurs résultats inhabituels.

2.5.2 Diagnostic par Baujat

Une méthode complémentaire a éré proposéepar Baujat et al: Voir papier Baujat B, Mahé C, Pignon JP, Hill C (2002): A graphical method for exploring heterogeneity in meta-analyses: Application to a meta-analysis of 65 trials. Statistics in Medicine, 30, 2641–52

baujat(ma_model_5)

-> Interprétation des sorties

L’axe X représente la contribution de l’essai au test Cochran Q global pour l’hétérogénéité
L’axe des Y représente l’influence de l’essai, définie comme la différence quadratique normalisée entre les effets du traitement estimés avec et sans l’essai.

2.5.3 Leave-one-out

on va retirer une à une les études présentes dans la base de données et re-calculer le summary effect size. On va pouvoir ainsi déterminer si certaines des études ont une large influence sur le modèle global.

LOO<-data.frame(leave1out(ma_model_4))
LOO$study<-row.names(LOO)
ggplot(LOO)+ geom_point(aes(study,estimate))+
  geom_errorbar(aes(study,ymin=ci.lb ,ymax=ci.ub ))+
  coord_flip()+
  theme_bw()

-> Interprétation des sorties

On remarque que 2 etudes impactent très fortement les résultats de la méta-analyse. Etudes à vérifier précisément.

2.6 Analysis of the publication bais

Un biais de publication désigne en science le fait que les chercheurs et les revues scientifiques ont bien plus tendance à publier des expériences ayant obtenu un résultat positif (statistiquement significatif) que des expériences ayant obtenu un résultat négatif (soutenant l’hypothèse nulle). Ce biais de publication donne aux lecteurs une perception biaisée (vers le positif) de l’état de la recherche Il est possible d’identifier dans un ensemble d’études, s’il y a présence ou pas de biais de publication

La présence d’un biais de publication (ou plus exactement, l’asymétrie de la courbe en entonnoir (funnel plot) ou les “effets de petites études”) et son impact potentiel sur les résultats peuvent être examinés par diverses méthodes, notamment :

1. le test de corrélation de rang (fonction ranktest()),
2. Test de régression d’Egger (fonction regtest()),
3. la méthode trim and fill (fonction trimfill()),
4. l’approche Henmi et Copas (fonction hc()),
5. D’autres méthodes comme les méthodes Rosenthal, Orwin et Rosenberg (fonction fsn()).

Nous allons focaliser dans ce TD sur deux méthodes.

2.6.1 Funnel plot

Il suppose que les petites études sont plus susceptibles d’être sujettes à un biais de publication que les grandes, et c’est cette différence qui est détectable. Si un chercheur termine un essai randomisé de grande envergure, il est probable qu’il voudra le voir publié même si le résultat est négatif en raison de l’effort requis. Pour les petits procès, cependant, la situation peut être différente. S’il existe un biais de publication, il est fort probable qu’il soit dû au fait que de petits essais négatifs n’ont pas été publiés.

Le diagramme en entonnoir est une représentation graphique de la taille des essais par rapport à l’ampleur de l’effet qu’ils rapportent. Comme la taille de l’essai augmente, les essais sont susceptibles de converger vers la taille réelle de l’effet sous-jacent. On s’attendrait à voir une dispersion égale des essais de part et d’autre de cet effet sous-jacent réel. Lorsqu’il y a eu biais de publication, on s’attend à une asymétrie dans la dispersion des petites études, davantage d’études ayant donné un résultat positif que celles ayant donné un résultat négatif.

funnel(ma_model_1)

-> interprétation des résultats

Il ne semble pas y avoir de biais depublication dans les études de ce jeu de données

2.6.2 trim-and-fill analysis

res.tf <- trimfill(ma_model_4)
funnel(res.tf, legend=TRUE, cex=1.2)

-> interprétation des résultats

Etant donné l’absence de biais de publication, cette analyse n’aboutit pas.