Metanálisis en R

Martín Alonso Rondón Sepúlveda http://rpubs.com/mrondon/ (Pontificia Universidad Javeriana)https://www.javeriana.edu.co

Metanálisis para Datos Binarios en R

Bases de datos:

Para este ejemplo, utilizaremos una revisión sistemática de 11 ensayos clínicos que compararon el Haloperidol versus placebo en el tratamiento de los síntomas en la esquizofrenia (Adams & Lawrie, 2013). El resultado clínico de interés es mejoría, por lo tanto, Riesgos Relativos (RR) mayores a 1 indican que el Haloperidol es mejor en comparación con el placebo. Los siguientes son los datos:

base <- data.frame(
  Estudio = c("Bechelli", "Garcia", "Kane", "Durost", "Garry", "Howard", "Selman", 
              "Serafetinides", "Simpson", "Spencer","Vichaiya"),
  año = c(1983, 2009, 2010, 1964, 1962, 1974, 1976, 1972, 1967, 1992, 1971),
  exitoshalo = c(9, 17, 63, 14, 19, 9, 9, 10, 6, 3, 20),
  fracasoshalo = c(21, 41, 49, 5, 7, 8, 20, 4, 10, 9, 9),
  totalhalo = c(30, 58, 112, 19, 26, 17, 29, 14, 16, 12, 29),
  exitosplac = c(25, 36, 81, 15, 22, 10, 16, 14, 5, 11,29),
  fracasosplac = c(6, 25, 41, 1, 4, 3, 13, 1, 2, 1,1),
  totalplac = c(31, 61, 122, 15, 26, 13, 29, 14, 7, 12,30),
  tiemposem = c("<6", "<6", "<6", "≥6", "≥6", "≥6", "≥6", "≥6", "≥6", "≥6", "≥6")
)

Estudio: autor del estudio.

año: año de publicación

exitoshalo: Número de participantes que mejoraron en el grupo de haloperidol

fracasoshalo: Número de participantes que no mejoraron en el grupo de haloperidol

totalhalo: Número de participantes en el grupo de haloperidol

exitosplac: Número de participantes que mejoraron en el grupo placebo.

fracasosplac: Número de participantes que no mejoraron en el grupo placebo

totalplac: Número de participantes en el grupo placebo

tiemposem: Tiempo de seguimiento en el estudio (en semanas)

El siguiente comando se utiliza para fijar la base para los análisis:

attach(base)

Metanálisis de efectos fijos (Mantel-Haenszel), reportando Riesgos Relativos (RR)

library(meta) 
library(metafor)
base <- base[order(base$año),] # ordenando la base de datos por año de publicación
metaEF <- metabin(exitoshalo, totalhalo, exitosplac, totalplac, 
                  data=base,  sm="RR", method="MH", studlab=paste(Estudio, año),  comb.fixed = TRUE,
                  comb.random = FALSE)  
print(metaEF, digits = 4, only.fixed = TRUE)
Number of studies: k = 11
Number of observations: o = 722 (o.e = 362, o.c = 360)
Number of events: e = 443

                        RR           95%-CI     z  p-value
Common effect model 0.6740 [0.5992; 0.7582] -6.57 < 0.0001

Quantifying heterogeneity (with 95%-CIs):
 tau^2 = 0.0160 [0.0000; 0.2732]; tau = 0.1266 [0.0000; 0.5226]
 I^2 = 39.8% [0.0%; 70.3%]; H = 1.29 [1.00; 1.84]

Test of heterogeneity:
     Q d.f. p-value
 16.60   10  0.0837

Details of meta-analysis methods:
- Mantel-Haenszel method
- Restricted maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Calculation of I^2 based on Q
- Continuity correction of 0.5 in studies with
  zero cell frequencies

Haciendo el Forest Plot

forest(metaEF)

Metanálisis de efectos aleatorios (Mantel-Haenszel), reportando Riesgos Relativos (RR)

library(meta) 
library(metafor)
base <- base[order(base$año),] # ordenando la base de datos por año de publicación
metaEA <- metabin(exitoshalo, totalhalo, exitosplac, totalplac, 
                  data=base,  sm="RR", method="MH", studlab=paste(Estudio, año),  comb.fixed =F)
print(metaEA)
Number of studies: k = 11
Number of observations: o = 722 (o.e = 362, o.c = 360)
Number of events: e = 443

                         RR           95%-CI     z  p-value
Random effects model 0.6917 [0.6023; 0.7943] -5.22 < 0.0001

Quantifying heterogeneity (with 95%-CIs):
 tau^2 = 0.0160 [0.0000; 0.2732]; tau = 0.1266 [0.0000; 0.5226]
 I^2 = 39.8% [0.0%; 70.3%]; H = 1.29 [1.00; 1.84]

Test of heterogeneity:
     Q d.f. p-value
 16.60   10  0.0837

Details of meta-analysis methods:
- Inverse variance method
- Restricted maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Calculation of I^2 based on Q
- Continuity correction of 0.5 in studies with
  zero cell frequencies

Haciendo el Forest Plot

forest(metaEA)

Haciendo el Gráfico de Embudo

funnel(metaEF, comb.fixed =F)

Para el análisis de subgrupos según el tiempo de seguimiento se debe hacer el siguiente comando

metaEF <- metabin(exitoshalo, totalhalo, exitosplac, totalplac, data=base,  
                  sm="RR", method="I", studlab=paste(Estudio, año), comb.fixed =F)

metaEF2 <- update(metaEF, byvar=ifelse(tiemposem=="<6","<6", ">=6"),
                  print.byvar=FALSE)

Haciendo el Forest Plot

forest(metaEF2)

Para el análisis de subgrupos según el año de publicación (antes del 1990 y de 1990 en adelante) se debe hacer el siguiente comando

metaEF3 <- update(metaEF, byvar=ifelse(año<1990,"<1990", ">=1990"), 
                  print.byvar=FALSE, print.CMH=TRUE)

forest(metaEF3)

Gráfico de L’Abbe

Cada uno de los puntos en el gráfico representa el riesgo relativo correspondiente a los diferentes estudios, de modo que la diagonal que divide el gráfico en dos secciones dejando a uno de los lados los estudios favorables al grupo de tratamiento y al otro los favorables al grupo control. La presencia de puntos dispersos, que no se sitúen de forma paralela a dicha diagonal, indicará posible heterogeneidad.

labbe(metaEF3, xlab="Tasa de eventos (Placebo)", ylab="Tasa de 
eventos (Haloperidol)", studlab=TRUE, cex.lab= 1.3, cex.axis=1.2, 
main="Haloperidol vs Placebo en el tratamiento \n 
de los síntomas en la esquizofrenia", cex.main=1.4, bg="green") 

abline(a=0,b=1) 

La mayoría de los puntos se agrupan cerca de la diagonal o por debajo de ella,lo que indica que la tasa de eventos con Haloperidol fue menor en comparación con el placebo en todos los estudios. La dispersión de los puntos en torno a la línea punteada sugiere que hay variabilidad en la magnitud del efecto entre los estudios, ya que en algunos, la reducción es grande, mientras que en otros es menor.

Gráfico de Galbraith

Para un modelo de efectos fijos, el gráfico muestra el inverso de los errores estándar en el eje horizontal contra los tamaños de efecto observados o los resultados estandarizados por sus errores estándar correspondientes en el eje vertical.

library(meta)
galbraith(metaEF3, level = 0.95, col = "blue", lty = 2, 
          lwd = 2, xlab = "Precisión", ylab = "Efecto estandarizado")

Sin embargo, los gráficos de Galbraith y de L’Abbé sugieren cierto grado de heterogeneidad, dado que los puntos no están completamente simétricos alrededor de la línea central. Además, los puntos parecen estar más concentrados en la parte inferior izquierda, lo que sugiere una posible asimetría. Adicionalmente, el gráfico de Galbraith muestra que uno de los estudios cae fuera de las bandas de confianza en el primero (aquel que proporciona una menor estimación del efecto).

Sesgos de publicación

Además del grafico de embudo, exiten otras alternativas para evaluar el sesgo de publicación, tales como:

Prueba de Egger: Evalúa la asimetría del gráfico de embudo. Prueba de Begg: También evalúa el sesgo de publicación, pero es menos común y menos potente que la prueba de Egger. Prueba de regresión de Harbord: Una alternativa a la prueba de Egger en el caso de razones de riesgos (RR) y odds ratios (OR). Gráficos de Cortina y Gráficos de Baujat: Permiten evaluar la influencia de cada estudio en la estimación global del efecto.

Prueba de Egger

Representar la recta de regresión entre la precisión de los estudios (variable independiente) y el efecto estandarizado (variable dependiente). Cuando no hay sesgo de publicación la recta de regresión se origina en el cero del eje Y. Cuánto más se aleje del cero, mayor evidencia de sesgo de publicación.

metabias(metaEF, plotit=T)

Linear regression test of funnel plot asymmetry

Test result: t = -4.48, df = 9, p-value = 0.0015
Bias estimate: -2.4322 (SE = 0.5433)

Details:
- multiplicative residual heterogeneity variance (tau^2 = 0.5716)
- predictor: standard error
- weight:    inverse variance
- reference: Egger et al. (1997), BMJ

Visualmente, no se aptrecia una indicación fuerte de sesgo de publicación, sin embargo, al evaluar el valor de p encontramos que si hay diferencias estadísticamente significativas de asimetría en el gráfico de embudo (p=0.0015). La estimación del sesgo de -2.4322 respalda esta conclusión, sugiriendo que los estudios más pequeños podrían estar sobreestimando los efectos observados.

Prueba de Begg

Evalúa la presencia de asociación entre las estimaciones de los efectos y sus varianzas. En caso de existir correlación es porque el efecto de la intervención es dependiente de la heterogeneidad. No es confiable tampoco con pocos estudios.

metabias(metaEF, method.bias="rank", plotit=T)

Rank correlation test of funnel plot asymmetry

Test result: z = -2.72, p-value = 0.0064
Bias estimate: -35.0000 (SE = 12.8452)

Reference: Begg & Mazumdar (1993), Biometrics

Se observa que los estudios con mayor varianza tienden a tener efectos extremos presentando una mayor dispersión en sus efectos estimados. Igual que en el grafico anterior, vemos que hay evidencia significativa de asimetría (p=0.0064), indicando posible sesgo de publicación.

Prueba de Harbord

Realiza la prueba de Harbord, que es adecuada para metaanálisis de razones de riesgos o odds ratios.

metabias(metaEF, method.bias = "harbord", plotit=T)

Linear regression test of funnel plot asymmetry

Test result: t = -3.28, df = 9, p-value = 0.0096
Bias estimate: -2.5368 (SE = 0.7742)

Details:
- multiplicative residual heterogeneity variance (tau^2 = 1.0049)
- predictor: standard error of score
- weight:    inverse variance of score
- reference: Harbord et al. (2006), Stat Med

Nuevamente encontramos que hay evidencia significativa de asimetría en el gráfico de embudo, lo que sugiere la posible presencia de sesgo de publicación. Nuevamente una estimación del sesgo negativa (-2.4163) indica que los estudios más pequeños podrían estar sobreestimando los efectos observados.

Gráfico de Cortina

Permite evaluar la influencia de cada estudio en la estimación global del efecto en un metanálisis, representando el valor p en función del estimador del efecto para cada estudio. Cada curva muestra cómo varía el valor p según el intervalo de confianza, mientras que líneas punteadas indican niveles de significancia (p = 0.01, 0.05, 0.10) y los intervalos de confianza (90%, 95%, 99%). Este gráfico es útil para evaluar la robustez del metanálisis, identificar estudios con alta influencia en los resultados y detectar posibles sesgos de publicación o heterogeneidad. Si las curvas están bien centradas y dentro de los intervalos esperados, se puede concluir que el metanálisis es sólido.

drapery(metaEF, 
        labels = "studlab",
        type = "pval", 
        legend = FALSE)

El gráfico de cortina muestra la relación entre el Risk Ratio y el valor p para cada estudio incluido en el metanálisis. La mayoría de los estudios presentan estimaciones agrupadas cercanas a 1.0, lo que sugiere un efecto combinado estable y consistente. La curva azul, que representa la estimación global, es estrecha, lo que indica una buena precisión en la estimación del efecto. Además, la mayoría de los estudios se encuentran dentro de los intervalos de confianza del 90%, 95% y 99%, lo que refuerza la solidez del metanálisis.

Sin embargo, es importante ver que algunos estudios muestran valores p más altos y curvas que se desvían de la tendencia general, lo que podría sugerir cierta heterogeneidad o la presencia de estudios con menor influencia en la estimación final.

Diagnóstico de influencia

Método de Baujat

Identifica estudios con alta influencia y heterogeneidad en el meta-análisis, mostrando una relación entre la contribución de cada estudio a la heterogeneidad y su impacto en el resultado. Útil para detectar estudios atípicos que generan heterogeneidad.

if (!require("remotes")) {
  install.packages("remotes")
}
remotes::install_github("MathiasHarrer/dmetar")
No outliers detected (random-effects model).
m.gen.inf <- InfluenceAnalysis(metaEF, random = TRUE)
[===========================================================================] DONE 
plot(m.gen.inf, "baujat")

Del anterior análisis, se puede concluir que los estudios Kane 2010, Bechelli 1983 y Garcia 2009 son los que más influyen en la heterogeneidad del metanálisis, por lo que se recomienda realizar un análisis de sensibilidad excluyendo estos estudios para evaluar su impacto en los resultados.

Datos Atípicos e influyentes

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

El gráfico de influencia muestra que los estudios Kane 2010, Bechelli 1983 y Garcia 2009 son los que más influyen en la estimación global del efecto, ya que se encuentran más alejados de la línea de referencia. Esto sugiere que estos estudios podrían estar sesgando los resultados del metanálisis, por lo que se recomienda realizar un análisis de sensibilidad excluyendo estos estudios para evaluar su impacto en los resultados.

DIAGRAMA DE ANÁLISIS DE SENSIBILIDAD

Determina la influencia de cada uno de los estudios en la estimación global del efecto. Consiste en realizar el metanálisis tantas veces como estudios seleccionados, de forma que cada vez se omite un estudio. Si los resultados de los distintos metanálisis son similares, se puede concluir que el resultado del metanálisis tiene solidez.

Análisis de influencia.

Las estimaciones agrupadas se calculan omitiendo un estudio a la vez.

meta1<-metainf(metaEF3, pooled="random") 

forest(meta1, text.random="Modelo de efectos aleatorios", 
comb.random=TRUE, leftlabs=c("año"), pooled.totals=F, 
pooled.events=F, xlab="Riesgo Relativo (RR)", fs.xlab=14,ff.xlab=2, 
col.square="red", col.diamond="blue") 

Al omitir estudios, los valores de RR varían entre 0.66 y 0.73, lo que sugiere que ningún estudio tiene un impacto desproporcionado en los resultados generales. Además, los valores de Tau² y Tau muestran leves variaciones, lo que indica que la heterogeneidad no cambia significativamente al excluir estudios individuales.

Los valores de P-value son menores a 0.0001, lo que indica que los resultados fueron estadísticamente significativos en todos los casos, por lo tanto, podemos decir que el metanálisis es robusto, ya que la omisión de estudios no altera sustancialmente las conclusiones.

Metanálisis acumulativo.

Se añaden los estudios a medida que fueron apareciendo. Este gráfico permite evaluar la solidez en la información y, si se llegó a un equilibrio en el resultado.

metaEF4 <- metabin(exitoshalo, totalhalo, exitosplac, totalplac, data=base, 
                   sm="RR", method="I", studlab=paste(Estudio, año), print.CMH=TRUE)

meta2<-metacum(metaEF4, pooled="random") 

forest(meta2, xlab="Riesgo Relativo (RR)", text.random ="Modelo de efectos aleatorios", 
            fs.xlab=14, ff.xlab=2, col.square="red", col.diamond= "blue", leftlabs=c("año")) 

El gráfico muestra que a medida que se incluyen más estudios, el RR se estabiliza y se vuelve más preciso, alcanzando 0.69 [0.60; 0.79] con el conjunto completo de estudios. A partir del año 1992, los sucesivos estudios realizados no añadieron ninguna información, ya que el resultado del metanálisis fue el mismo.

Ejemplo con desenlaces continuos

Usando la base SuicidePrevention de la librería dmetar (Harrer et al., 2021), se realizará un metanálisis para desenlaces continuos.

if (!require("remotes")) {
  install.packages("remotes")
}
remotes::install_github("MathiasHarrer/dmetar")
library(dmetar)
library(esc)
library(tidyverse)
data(SuicidePrevention)
data(SuicidePrevention)
base <- SuicidePrevention
paged_table(base)
library(meta)
library(metafor)
mEAcont <- metacont(n.e = n.e, mean.e = mean.e, sd.e = sd.e,
                   n.c = n.c, mean.c = mean.c, sd.c = sd.c,
                   studlab = author, data = SuicidePrevention,
                   sm = "SMD", method.smd = "Hedges",
                   fixed = FALSE, random = TRUE,
                   method.tau = "REML", method.random.ci = "HK",
                   title = "Suicide Prevention")
summary(mEAcont)
Review:     Suicide Prevention

                    SMD             95%-CI %W(random)
Berry et al.    -0.1428 [-0.4315;  0.1459]       15.6
DeVries et al.  -0.6077 [-0.9402; -0.2752]       12.3
Fleming et al.  -0.1112 [-0.6177;  0.3953]        5.7
Hunt & Burke    -0.1270 [-0.4725;  0.2185]       11.5
McCarthy et al. -0.3925 [-0.7884;  0.0034]        9.0
Meijer et al.   -0.2676 [-0.5331; -0.0021]       17.9
Rivera et al.    0.0124 [-0.3454;  0.3703]       10.8
Watkins et al.  -0.2448 [-0.6848;  0.1952]        7.4
Zaytsev et al.  -0.1265 [-0.5062;  0.2533]        9.7

Number of studies: k = 9
Number of observations: o = 1147 (o.e = 571, o.c = 576)

                         SMD             95%-CI     t p-value
Random effects model -0.2304 [-0.3734; -0.0874] -3.71  0.0059

Quantifying heterogeneity (with 95%-CIs):
 tau^2 = 0.0044 [0.0000; 0.0924]; tau = 0.0661 [0.0000; 0.3040]
 I^2 = 7.4% [0.0%; 67.4%]; H = 1.04 [1.00; 1.75]

Test of heterogeneity:
    Q d.f. p-value
 8.64    8  0.3738

Details of meta-analysis methods:
- Inverse variance method
- Restricted maximum-likelihood estimator for tau^2
- Q-Profile method for confidence interval of tau^2 and tau
- Calculation of I^2 based on Q
- Hartung-Knapp adjustment for random effects model (df = 8)
- Hedges' g (bias corrected standardised mean difference;
  using exact formulae)

Haciendo el Forest Plot

forest(mEAcont)

Haciendo el Gráfico de Embudo

funnel(mEAcont, comb.fixed =F)

Para el análisis de subgrupos según el tiempo de seguimiento se debe hacer el siguiente comando

mEAcont2 <- update(mEAcont, byvar=control, print.byvar=FALSE)

Haciendo el Forest Plot

forest(mEAcont2)

DIAGRAMA DE ANÁLISIS DE SENSIBILIDAD

Determina la influencia de cada uno de los estudios en la estimación global del efecto. Consiste en realizar el metanálisis tantas veces como estudios seleccionados, de forma que cada vez se omite un estudio. Si los resultados de los distintos metanálisis son similares, se puede concluir que el resultado del metanálisis tiene solidez.

Análisis de influencia.

Las estimaciones agrupadas se calculan omitiendo un estudio a la vez.

meta1<-metainf(mEAcont2, pooled="random") 

forest(meta1, text.random="Modelo de efectos aleatorios", 
comb.random=TRUE, leftlabs=c("año"), pooled.totals=F, 
pooled.events=F, xlab="Riesgo Relativo (RR)", fs.xlab=14,ff.xlab=2, 
col.square="red", col.diamond="blue") 

Metanálisis acumulativo.

Se añaden los estudios a medida que fueron apareciendo. Este gráfico permite evaluar la solidez en la información y, si se llegó a un equilibrio en el resultado.

meta2<-metacum(mEAcont2, pooled="random") 

forest(meta2, xlab="Riesgo Relativo (RR)", text.random ="Modelo de efectos aleatorios", 
            fs.xlab=14, ff.xlab=2, col.square="red", col.diamond= "blue", leftlabs=c("año")) 

El gráfico muestra que a partir del año 1992, los sucesivos estudios realizados no añadieron ninguna información, ya que el resultado del metanálisis fue el mismo.

Adams, B., CE, & Lawrie, S. (2013). Haloperidol versus placebo for schizophrenia. Cochrane Database of Systematic Reviews, 11. https://doi.org/10.1002/14651858.CD003082.pub3
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.

References