Untitled

Comenzamos cargandos las librerias necesarias que son metafor y el tidyverse, instalandolas previamente si aún no lo están:

if (!require("metafor")) install.packages("metafor")

Loading required package: metafor

Loading required package: Matrix

Loading required package: metadat

Loading required package: numDeriv

Loading the 'metafor' package (version 4.8-0). For an
introduction to the package please type: help(metafor)

if (!require("tidyverse")) install.packages("tidyverse")

Loading required package: tidyverse

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand() masks Matrix::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ tidyr::pack()   masks Matrix::pack()
✖ tidyr::unpack() masks Matrix::unpack()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(metafor)
library(tidyverse)

En primer lugar cargamos una hoja del fichero excel con las frecuencias de interes en el formato que queremos:

hojaActual="OKFET"
datos <- readxl::read_excel("Cualitativas_OK.xlsx", sheet = hojaActual)
datos

# A tibble: 15 × 18
   Autor     N N_hombre N_ictus N_SCIPara N_ventilacionprolong N_isquemInt N_DRA
   <chr> <dbl>    <dbl>   <dbl>     <dbl>                <dbl>       <dbl> <dbl>
 1 Arno…    74       29       6         2                   14          NA    NA
 2 Arno…    49       26       4         2                   13          NA    NA
 3 Di E…    21       18       2         1                    4           0    NA
 4 Kim …    55       42       2         1                    1           1    NA
 5 Kim …    20       15       2         0                    0           1    NA
 6 Koiz…    60       51       2         2                    2          NA     2
 7 Liak…    36       21       2         0                   14           1     8
 8 Panf…    26       22       0         1                    5          NA    NA
 9 Scho…    32       18       2         4                    3          NA    NA
10 Scho…    68       47       2         5                    2          NA    NA
11 Shre…    82       NA       2         5                   14          NA    16
12 Shre…    69       NA       2         3                   18          NA    15
13 Tian…   698      530       2        35                   NA          NA   102
14 Verh…    94       62       2         4                   NA          NA    23
15 Weis…    41       32       4         2                   10          NA    NA
# ℹ 10 more variables: N_dialisis <dbl>, N_dialisistrans <dbl>,
#   N_dialisisperm <dbl>, N_sangrado_retoracotomía <dbl>,
#   N_reintervencion <dbl>, N_paralisisrecurrente <dbl>, N_endoleak <dbl>,
#   N_mortintrahos <dbl>, N_mort30días <dbl>, N_mortconjunta <dbl>

Realmente los conteos que tenemos disponibles son estos:

 datos %>% select(starts_with("N_")) %>% names() %>% paste0(sep=" ", collapse=" ") %>% cat()

N_hombre  N_ictus  N_SCIPara  N_ventilacionprolong  N_isquemInt  N_DRA  N_dialisis  N_dialisistrans  N_dialisisperm  N_sangrado_retoracotomía  N_reintervencion  N_paralisisrecurrente  N_endoleak  N_mortintrahos  N_mort30días  N_mortconjunta 

Generamos un informe para la columna:

columnaInteres="N_hombre"
Descripcion="Hombres"

datos$Ncasos=datos[[columnaInteres]]

datosCalculo=datos %>% select(Autor, N, Ncasos) %>% filter(complete.cases(.))
nEstudios=nrow(datosCalculo)

La forma de trabajar es fijarnos en la columna Autor, N y Ncasos (columna de interés). Queremos hacer una metaanálisis de porcentajes calculados como Ncasos/N. Para ello usamos la función escalc del paquete metafor, que calcula los efectos y sus varianzas. En este caso usamos la transformación logit para que los intervalos de confianza estén acotados entre 0 y 1.

datosAmpliados <- escalc(measure="PLO", xi=Ncasos, ni=N, data=datosCalculo, append=TRUE)
summary(datosAmpliados) %>% as_tibble()

# A tibble: 13 × 10
   Autor           N Ncasos     yi      vi    sei     zi     pval   ci.lb  ci.ub
   <chr>       <dbl>  <dbl>  <dbl>   <dbl>  <dbl>  <dbl>    <dbl>   <dbl>  <dbl>
 1 Arnold et …    74     29 -0.439 0.0567  0.238  -1.85  6.50e- 2 -0.906  0.0274
 2 Arnold et …    49     26  0.123 0.0819  0.286   0.428 6.68e- 1 -0.438  0.684 
 3 Di Eusanio…    21     18  1.79  0.389   0.624   2.87  4.06e- 3  0.570  3.01  
 4 Kim et al …    55     42  1.17  0.101   0.317   3.69  2.20e- 4  0.551  1.79  
 5 Kim et al …    20     15  1.10  0.267   0.516   2.13  3.34e- 2  0.0865 2.11  
 6 Koizumi et…    60     51  1.73  0.131   0.362   4.80  1.61e- 6  1.03   2.44  
 7 Liakopoulo…    36     21  0.336 0.114   0.338   0.995 3.20e- 1 -0.326  0.999 
 8 Panfilov e…    26     22  1.70  0.295   0.544   3.14  1.71e- 3  0.639  2.77  
 9 Schoeberl …    32     18  0.251 0.127   0.356   0.705 4.81e- 1 -0.447  0.950 
10 Schoeberl …    68     47  0.806 0.0689  0.262   3.07  2.15e- 3  0.291  1.32  
11 Tian et al    698    530  1.15  0.00784 0.0885 13.0   1.67e-38  0.975  1.32  
12 Verhoye et…    94     62  0.661 0.0474  0.218   3.04  2.38e- 3  0.235  1.09  
13 Weiss et al    41     32  1.27  0.142   0.377   3.36  7.74e- 4  0.529  2.01  

Ya tenemos los datos preparados para hacer el metaanálisis. Usamos la función rma para hacer el metaanálisis, en este caso usando un modelo de efectos aleatorios (REML es el método por defecto).

res <- rma(yi, vi, data=datosAmpliados, method="REML")
summary(res)

Random-Effects Model (k = 13; tau^2 estimator: REML)

  logLik  deviance       AIC       BIC      AICc   
-12.4421   24.8843   28.8843   29.8541   30.2176   

tau^2 (estimated amount of total heterogeneity): 0.3365 (SE = 0.1842)
tau (square root of estimated tau^2 value):      0.5801
I^2 (total heterogeneity / total variability):   82.76%
H^2 (total variability / sampling variability):  5.80

Test for Heterogeneity:
Q(df = 12) = 65.7433, p-val < .0001

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub      
  0.8284  0.1876  4.4163  <.0001  0.4608  1.1961  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ahora mostramos los resultados en un forest plot.

Vamos a posicionar la columna de casos algo más a la izquierda para que no se solape con los intervalos de confianza. Para ello usamos el argumento ilab.xpos. Podemos añadir más columnas si queremos, por ejemplo la columna N:

forest(res, slab=datosCalculo$Autor, xlab="Proporción", atransf=plogis, ilab=cbind(datosCalculo$Ncasos, datosCalculo$N), ilab.xpos=c(-3.5, -3))
#Poner los titulos de las columnas
text(-3.5, nEstudios+0.75, "Casos")
text(-3, nEstudios+0.75, "N")

op <- par(cex=0.75, font=2)

AHora podríamos estudiar el sesgo de publicación con un funnel plot y el test de Egger:

funnel(res, atransf=plogis)

regtest(res, model="rma", predictor="sei")

Regression Test for Funnel Plot Asymmetry

Model:     mixed-effects meta-regression model
Predictor: standard error

Test for Funnel Plot Asymmetry: z = 1.4543, p = 0.1459
Limit Estimate (as sei -> 0):   b = 0.1664 (CI: -0.7974, 1.1301)