Application 1 : NSAIDS and acute pain

Exemple tiré de (Schwarzer, Carpenter, & Rücker, 2015, p. 108).

Moore, Tramer, Carroll, Wiffen, & McQuay (1998) conducted a systematic review of 37 randomized placebo-controlled trials on the effectiveness and safety of topical non-steroidal anti-inflammatory drugs (NSAIDS) in acute pain. The main outcome was treatment success, defined as a reduction in pain of at least 50 %.

Les données sont disponibles dans le fichier dataset11.csv.

d <- read.csv("data/dataset11.csv")
d

1. À l’aide d’une analyse FE couplée à une analyse RE sur la base des odds-ratio calculés à partir des résultats disponibles dans le fichier, retrouver les résultats indiqués dans le graphique ci-dessus. Isoler les études présentant les effets les plus extrêmes.
2. Réaliser un graphique de type “funnel plot” à l’aide de la commande funnel. Quelle conclusion peut-on tirer de ces résultats ?
3. Utiliser un test non paramétrique de Begg and Mazumdar et conclure quant à l’hypothèse nulle de symétrie du funnel plot, en d’autres termes de l’absence de bias de publication.

Application 2

Exemple adapté de (Schwarzer et al., 2015, p. 152).

Il s’agit des données exposées dans l’application n° 1 du premier TP. On rappelle qu’il s’agit d’une étude sur les bronchodilatateurs et que la variable réponse est l’évolution du volume expiratoire de pointe au cours du suivi. Seules les 7 études présentant les erreurs standard les plus faibles ont été retenues pour cette analyse, et des données manquantes ont été artificiellement générées pour certains des centres et les nouveaux effectifs observés sont désignés par les variables Neo et Nco (cf. variables Nem et Ncm qui indiquent le nombre total de valeurs manquantes présentes dans chaque étude).

Les données sont disponibles dans le fichier dataset01m.rda.

load("data/dataset01m.rda")
d

##                 study Ne    Me    Se Nc    Mc    Sc Nem Ncm Neo Nco
## 2          Boner 1989 20 15.70 13.10 20 22.70 16.47   0   0  20  20
## 5  DeBenedictis 1994a 17 14.40 11.10 17 27.40 17.30   0   0  17  17
## 12     Novembre 1994f 24 15.42  8.35 24 28.46 13.84   8   8  16  16
## 13     Novembre 1994s 19 11.00 12.40 19 26.10 14.90   0   0  19  19
## 14         Oseid 1995 20 14.10  9.50 20 28.90 18.00   5   5  15  15
## 16          Shaw 1985  8 10.27  7.02  8 34.43 10.96   0   0   8   8
## 17        Todaro 1993 13 10.10  8.90 13 23.50  4.00   0   0  13  13

Dans un premier temps, voici les résultats d’une méta-analyse sur les 7 études, en considérant les effectifs modifiés :

m <- metacont(Neo, Me, Se, Nco, Mc, Sc, data = d, studlab = study)
print(m, digits = 3)

##                         MD             95%-CI %W(fixed) %W(random)
## Boner 1989          -7.000 [-16.223;   2.223]      10.8       11.8
## DeBenedictis 1994a -13.000 [-22.771;  -3.229]       9.6       10.7
## Novembre 1994f     -13.040 [-20.960;  -5.120]      14.7       15.2
## Novembre 1994s     -15.100 [-23.816;  -6.384]      12.1       13.0
## Oseid 1995         -14.800 [-25.100;  -4.500]       8.7        9.7
## Shaw 1985          -24.160 [-33.179; -15.141]      11.3       12.2
## Todaro 1993        -13.400 [-18.704;  -8.096]      32.7       27.4
## 
## Number of studies combined: k = 7
## 
##                           MD             95%-CI     z  p-value
## Fixed effect model   -14.161 [-17.196; -11.127] -9.15 < 0.0001
## Random effects model -14.222 [-17.652; -10.792] -8.13 < 0.0001
## 
## Quantifying heterogeneity:
## tau^2 = 3.8428; H = 1.10 [1.00; 1.62]; I^2 = 17.9% [0.0%; 61.7%]
## 
## Test of heterogeneity:
##     Q d.f. p-value
##  7.31    6  0.2934
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2

Etant donné le contexte, on s’intéresse essentiellement au modèle à effet aléatoire. Les paramètres estimés peuvent être retrouvés à partir de la sortie ci-dessus ou directement à partir de l’objet m :

cat("TE=", m$TE.random, ", 95% CI: ", m$lower.random, ";", m$upper.random, "\n", sep = "")

## TE=-14.22221, 95% CI: -17.65205;-10.79238

Pour réaliser un diagramme de type “funnel plot”, on utilise la commande funnel. La syntaxe du modèle et de la fonction funnel sont sensiblement les mêmes avec le package metafor (ainsi que les sorties, à l’exception du fait que metafor::rma ne calcule qu’un modèle à la fois).

funnel(m, text = d$study, cex = 0.7)

Pour examiner l’impact des données manquantes on va procéder à l’aide d’une analyse en deux étapes (White, Higgins, & Wood, 2008) : ré-estimation des effets et de leur erreur standard après imputation puis estimation des paramètres du modèle à effet aléatoire. Un premier choix d’analyse de sensibilité pour les données manquantes est de considérer que $\mu_e = \mu_c = \mu$, avec $\nu_e^2 = \nu_c^2 = 0$ (stratégie dite “fixed equal”). En d’autres termes l’effet $\delta$ est le même entre les deux groupes pour chaque étude et entre les études. Cette analyse est sensible au déséquilibre des valeurs manquantes entre les groupes. Dans le cas présent, cela ne pose pas de problème car le nombre d’observations manquantes est équilibré entre les deux groupes. Un autre choix d’analyse de sensibilité consiste à postuler que les différences sont constantes entre les deux groupes, mais varient d’une étude à l’autre (stratégie dite “random equal”). Dans ce cas, on considère donc que $\delta_e = \delta_c = \delta$, avec $\delta \sim N(\mu,\nu^2)$. Les études présentant un déséquilibre en termes de valeurs manuantes se verront attribuées un poids plus faible avec cette stratégie.

La fonction auxiliaire ci-dessous permet de calculer l’erreur standard de l’estimateur $\hat\Delta$ qui reflète l’effet moyen classique (sur données complètes) auquel on incorpore un effet partiel ajusté sur la proportion de données manquantes dans les deux groupes (Schwarzer et al., 2015, p. 145) :

semiss <- function(se, n.e, p.e, n.c, p.c, mu.e, nu.e, mu.c, nu.c, rho) {
  V1 <- function(p.e, p.c, nu.e, nu.c, rho)
    nu.e^2*p.e^2 + nu.c^2*p.c^2 - 2*rho*nu.e*nu.c*p.e*p.c
  V2 <- function(n.e, p.e, n.c, p.c, mu.e, nu.e, mu.c, nu.c) 
    (mu.e^2+nu.e^2)*p.e*(1-p.e)/n.e + (mu.c^2+nu.c^2)*p.c*(1-p.c)/n.c
  sqrt(se^2 + V1(p.e, p.c, nu.e, nu.c, rho) + V2(n.e, p.e, n.e, p.c, mu.e, nu.e, mu.c, nu.c))
}

Voici en pratique ce que donnent ces deux stratégies sur les données :

## treatment effect from the initial model
d$TE <- m$TE 
d$seTE <- m$seTE
## fixed equal strategy: mu_e = mu_c = 8; nu_e = nu_c = 0
mu.e <- mu.c <- 8
nu.e <- nu.c <- 0
rho <- 0   ## not used but would reflect the within-study effect correlation
## proportion missing in each study arm
d$Pe <- with(d, Pe <- (Ne - Neo) / Ne)
d$Pc <- with(d, Pc <- (Nc - Nco) / Nc)
## mean effect for each study under strategy 1
d$TEs1 <- with(d, TE + mu.e*Pe - mu.c*Pc)
d$seTEs1 <- with(d, semiss(seTE, Neo, Pe, Nco, Pc, mu.e, nu.e, mu.c, nu.c, rho))
## dummy variable coding for studies with missings
selmiss <- d$study %in% c("Novembre 1994f", "Oseid 1995")

La procédure ci-dessus résulte en une légère augmentation des variances des estimateurs des effets (les moyennes restent inchangées car la proportion de valeurs manquantes dans les deux groupes est identique) :

d[selmiss, c("study", "TE", "seTE", "TEs1", "seTEs1")]

##             study     TE     seTE   TEs1   seTEs1
## 12 Novembre 1994f -13.04 4.040947 -13.04 4.255236
## 14     Oseid 1995 -14.80 5.255156 -14.80 5.405244

Et voici les résultats du modèle avec ces valeurs ajustées :

m1 <- metagen(TEs1, seTEs1, data = d, studlab = study, comb.fixed = FALSE)
cat("Fixed-equal TE=", m1$TE.random, ", 95% CI: ", m1$lower.random, ";", m1$upper.random, "\n", sep = "")

## Fixed-equal TE=-14.23488, 95% CI: -17.69984;-10.76993

Pour la seconde stratégie, voici la procédure correspondante :

mu.e <- mu.c <- 8
nu.e <- nu.c <- sqrt(5)
rho <- 1
d$TEs2 <- with(d, TE + mu.e*Pe - mu.c*Pc)
d$seTEs2 <- with(d, semiss(seTE, Neo, Pe, Nco, Pc, mu.e, nu.e, mu.c, nu.c, rho))
d[selmiss, c("study", "TE", "seTE", "TEs2", "seTEs2")]

##             study     TE     seTE   TEs2   seTEs2
## 12 Novembre 1994f -13.04 4.040947 -13.04 4.271525
## 14     Oseid 1995 -14.80 5.255156 -14.80 5.416795

Ainsi que les résultats du nouveau modèle de méta-régression :

m2 <- metagen(TEs2, seTEs2, data = d, studlab = study, comb.fixed = FALSE)
cat("Random-equal TE=", m2$TE.random, ", 95% CI: ", m2$lower.random, ";", m2$upper.random, "\n", sep = "")

## Random-equal TE=-14.23579, 95% CI: -17.70332;-10.76827

Les résultats sont similaires à ceux obtenus avec la première stratégie car les valeurs manquantes sont équilibrées entre les deux groupes. Le cas échéant, cette stratégie augmenterait encore plus l’erreur standard des études avec données manquantes en comparaison de la première stratégie.

Application 3 : Mucolytic agent and COPD

Exemple tiré de (Schwarzer et al., 2015, p. 41).

Poole & Black (2006) report a meta-analysis of mucolytic agents versus placebo for patients with chronic bronchitis and/or chronic obstructive pulmonary disease. The outcome is the mean number of acute exacerbations per month. Acute exacerbation is defined as an increase in cough and in the volume and/or purulence of sputum. As all studies report a mean number of exacerbations, we can work with mean differences, rather than standardised mean differences. R code to read in the data is given in Fig. 2.7. Notice that studies 5 and 12 (Jackson 1984, Grillage 1985) have zero standard errors.

Les données sont disponibles dans le fichier dataset03.csv.

d <- read.csv("data/dataset03.csv")
d

##         author year  Ne   Me   Se  Nc   Mc   Sc    duration
## 1   Bontognali 1991  30 0.70 3.76  30 1.27 4.58 <= 3 months
## 2  Castiglioni 1986 311 0.10 0.21 302 0.20 0.29 <= 3 months
## 3    Cremonini 1986  21 0.25 0.23  20 0.71 0.29 <= 3 months
## 4       Grassi 1994  42 0.16 0.29  41 0.45 0.43 <= 3 months
## 5      Jackson 1984  61 0.11 0.00  60 0.13 0.00 <= 3 months
## 6      Allegra 1996 223 0.07 0.11 218 0.11 0.14  > 3 months
## 7     Babolini 1980 254 0.13 0.18 241 0.33 0.27  > 3 months
## 8        Boman 1983  98 0.20 0.27 105 0.32 0.30  > 3 months
## 9       Borgia 1981  10 0.05 0.08   9 0.15 0.17  > 3 months
## 10    Decramer 2005 256 0.10 0.11 267 0.11 0.16  > 3 months
## 11      Grassi 1976  35 0.14 0.15  34 0.27 0.21  > 3 months
## 12    Grillage 1985  54 0.10 0.00  55 0.12 0.00  > 3 months
## 13      Hansen 1994  59 0.11 0.15  70 0.16 0.19  > 3 months
## 14     Malerba 2004 115 0.06 0.08 119 0.07 0.08  > 3 months
## 15     McGavin 1985  72 0.42 0.34  76 0.52 0.35  > 3 months
## 16     Meister 1986  90 0.15 0.15  91 0.20 0.19  > 3 months
## 17     Meister 1999 122 0.06 0.15 124 0.10 0.15  > 3 months
## 18     Moretti 2004  63 0.12 0.14  61 0.17 0.17  > 3 months
## 19       Nowak 1999 147 0.03 0.06 148 0.06 0.12  > 3 months
## 20    Olivieri 1987 110 0.18 0.31 104 0.33 0.41  > 3 months
## 21        Parr 1987 243 0.18 0.21 210 0.21 0.21  > 3 months
## 22        Pela 1999  83 0.17 0.18  80 0.29 0.32  > 3 months
## 23   Rasmussen 1988  44 0.13 0.21  47 0.14 0.19  > 3 months

1. Réaliser une analyse séparée (modèle FE et RE) dans le groupe défini par la variable duration. Comparer les effets estimés dans les deux sous-groupes.
2. Tester la différence entre les deux sous-groupes à l’aide de la fonction metagen (option sm = "MD") et conclure sur l’effet moyen dans les deux sous-groupes. Comparer avec les résultats produits par la commande metacont (option byvar = duration), en particulier en ce qui concerne le test d’hétérogénéité et le test de l’effet de groupe calculé avec metagen.