library(metafor)
#
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}

packages <- c("metafor", "esc", "readxl")
ipak(packages)

# http://www.metafor-project.org/doku.php

#
setwd("C:/Users/Desktop")
datosMETA <- read_excel("Meta-Analisis 17_6_22_zonas.xlsx", sheet="Hoja1", col_names = TRUE)
as.numeric(datosMETA$ALFA)
as.numeric(datosMETA$N)
datosMETA

res_ABT <- escalc(measure="ABT", ai=datosMETA$ALFA, mi=datosMETA$items, ni=datosMETA$N, data=datosMETA)
res_ABT

resultMA<-rma(yi=res_ABT, vi=res_ABT, 
              method="DL",
              weighted=TRUE, 
              level=95, 
              digits=4)

#
summary.rma(resultMA)
res_ABT
#
confint(resultMA)

#Media Total FS
ModeradorMediaTotalFS<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                           mods =~ res_ABT$Media.Total.FS, 
                           method="DL", 
                           weighted=TRUE, 
                           level=95, 
                           digits=4)
summary.rma(ModeradorMediaTotalFS)

#SD Total FS
ModeradorSDTotalFS<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                        mods =~ res_ABT$SD.Total.FS, 
                        method="DL", 
                        weighted=TRUE, 
                        level=95, 
                        digits=4)
summary.rma(ModeradorSDTotalFS)

#edad Media
Moderador.edadMedia<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                         mods =~ res_ABT$edad.Media, 
                         method="DL", 
                         weighted=TRUE, 
                         level=95, 
                         digits=4)
summary.rma(Moderador.edadMedia)

#edad SD
Moderador.Edad.SD<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                       mods =~ res_ABT$Edad.SD, 
                       method="DL", 
                       weighted=TRUE, 
                       level=95, 
                       digits=4)
summary.rma(Moderador.Edad.SD)

#N
Moderador.N<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                 mods =~ res_ABT$N, 
                 method="DL", 
                 weighted=TRUE, 
                 level=95, 
                 digits=4)
summary.rma(Moderador.N)

#N hombre
Moderador.N.homb<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                      mods =~ res_ABT$N.homb, 
                      method="DL", 
                      weighted=TRUE, 
                      level=95, 
                      digits=4)
summary.rma(Moderador.N.homb)

#N.muj
Moderador.N.muj<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                     mods =~ res_ABT$N.muj, 
                     method="DL", 
                     weighted=TRUE, 
                     level=95, 
                     digits=4)
summary.rma(Moderador.N.muj)

#muestra 
Moderador.muestra<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                       mods=~factor(res_ABT$muestra)-1, 
                       method="DL", 
                       weighted=TRUE, 
                       level=95, 
                       digits=4)
summary.rma(Moderador.muestra) 
table(datosMETA$muestra)

Moderador.muestra<- rma.uni(yi=res_ABT$yi, vi=res_ABT$vi,
                       mods = ~ factor(res_ABT$muestra),
                       method = "DL",
                       weighted = TRUE,
                       level = 95,
                       digits = 4)

summary.rma(Moderador.muestra)
#idioma
#
Moderador.IDIOMA<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                      mods=~factor(res_ABT$IDIOMA)-1, 
                      method="DL", 
                      weighted=TRUE, 
                      level=95, 
                      digits=4)
summary.rma(Moderador.IDIOMA) 
table(datosMETA$IDIOMA)

# Moderador categorico: país
#
Moderador.País<-rma(yi=res_ABT$yi, vi=res_ABT$vi, 
                    mods=~factor(res_ABT$ZONA.GEOGRAFICA)-1, 
                    method="DL", 
                    weighted=TRUE, 
                    level=95, 
                    digits=4)
summary.rma(Moderador.País) 

#
table(datosMETA$ZONA.GEOGRAFICA)

#
regtest(resultMA, model="lm")

#Rosenthal.
fsn(yi=res_ABT$yi, vi=res_ABT$vi,type = "Rosenthal")

#sensitividad:
leave1out(resultMA)
influenciado<-influence.rma.uni(resultMA)
influenciado$is.infl


##Forest plot
forest(resultMA,
       slab = datosMETA$Author,
       order = order(datosMETA$Author))

##Funnel Plot 
funnel(resultMA, level = c(90, 95, 99),
       shade = c("white", "gray55", "gray75"),
       refline = 0, legend = TRUE)

##Funnel Plot 
funnel(resultMA, level = c(90, 95, 99),
       shade = c("white", "gray55", "gray75"),
       legend = TRUE)