library(metafor)
library(meta)
library(tidyverse)


## Correct file extension: ".csv" instead of ".cvs"
##
prevalence_2020_cow_nomv  = read.csv("prev_cow.csv")


##
##
## 1) Random effects meta-analysis without subgroups
##
##

## Meta-analysis using metafor
##
ies.da = escalc(xi = nlameanimal, ni = ssizeanimal,
                data = prevalence_2020_cow_nomv,
                measure = "PFT", add = 0)
m4 = rma(yi, vi, data = ies.da, method = "DL")

## Meta-analysis using meta
##
m = metaprop(nlameanimal, ssizeanimal, author,
             data = ies.da,
             sm = "PFT", method.tau = "DL", method.ci = "NAsm",
             comb.fixed = FALSE, null.effect = 0)

## Exactly the same results for transformed prevalences
##
m4
predict(m4)
##
print(summary(m, backtransf = FALSE), digits.zval = 4)

## ... and for backtransformed proportions
##
p.m4 = predict(m4, transf = transf.ipft.hm,
               targ = list(ni = prevalence_2020_cow_nomv$ssizeanimal))
print(p.m4, digits = 4)
summary(m)





##
##
## 2) Random effects meta-analysis with subgroups - assuming common
##    between-study variance
##
##

## Subgroup meta-analysis with metafor
##
m4.lcmbi = rma(yi, vi, data = ies.da, method = "DL", mods = ~lcmbi)

## Subgroup meta-analysis with meta
##
m.lcmbi = update(m, byvar = lcmbi, tau.common = TRUE, print.byvar = FALSE)

## Again, exactly the same subgroup results for transformed prevalences
##
m4.lcmbi
print(summary(m.lcmbi, backtransf = FALSE), digits.zval = 4)
## Results for second subgroup (metafor)
pft2 = sum(coef(m4.lcmbi))
V = vcov(m4.lcmbi)
metagen(pft2, sqrt(sum(diag(V)) + 2 * V[1, 2]))

## Different results in meta and metafor
## (due to using a common harmonic mean in metafor and
##  subgroup-specific harmonic means in meta)
##
p.m4.lcmbi = predict(m4.lcmbi, transf = transf.ipft.hm,
                     targ = list(ni = prevalence_2020_cow_nomv$ssizeanimal))

print(p.m4.lcmbi, digits = 4)
summary(m.lcmbi)

## Overall harmonic mean (used in metafor)
##
n.harm = 1 / mean(1 / m.lcmbi$n)
n.harm
## Subgroup-specific harmonic means (used in meta)
##
n.harm.w = m.lcmbi$n.harmonic.mean.w
n.harm.w

## Prevalences in subgroups (metafor)
##
round(meta:::backtransf(m.lcmbi$TE.random.w, sm = "PFT", n = rep(n.harm, 2)), 4)

## Prevalences in subgroups (meta)
##
round(meta:::backtransf(m.lcmbi$TE.random.w, sm = "PFT", n = n.harm.w), 4)





##
##
## 3) Random effects meta-analysis with subgroups - allowing for
##    different between-study variances in subgroups
##
##

## Subgroup meta-analysis with metafor
##
m4.lcmbi.1 = rma(yi, vi, data = ies.da, method = "DL",
                 subset = lcmbi == "Locomotion Scoring Method")
m4.lcmbi.2 = rma(yi, vi, data = ies.da, method = "DL",
                 subset = lcmbi == "Records")
##
m4.lcmbi.1
m4.lcmbi.2

## Subgroup meta-analysis with meta
##
summary(update(m, byvar = lcmbi, tau.common = FALSE),
        print.byvar = FALSE, backtransf = FALSE)





##
##
## 4) Forest plot
##
##

pdf("metaprop04.pdf", width = 9, height = 10)
##
forest(m.lcmbi,
       xlim = c(0, 1),
       rightcols = FALSE,
       leftcols = c("studlab", "effect", "ci"),
       leftlabs = c("Study", "Prevalence", "95% C.I."),
       text.random = "Prevalence of Lameness in British Dairy Cattle",
       bylab = "Lameness Data Source",
       print.byvar = TRUE,
       xlab = "Prevalence of Lameness", smlab = "",
       squaresize = 0.5, col.square = "navy",
       col.diamond = "maroon", col.diamond.lines = "maroon",
       fs.hetstat = 10,
       print.Q = TRUE,
       print.pval.Q = TRUE,
       sortvar = seTE)
##
forest(update(m.lcmbi, tau.common = FALSE),
       xlim = c(0, 1),
       rightcols = FALSE,
       leftcols = c("studlab", "effect", "ci"),
       leftlabs = c("Study", "Prevalence", "95% C.I."),
       text.random = "Prevalence of Lameness in British Dairy Cattle",
       bylab = "Lameness Data Source",
       print.byvar = TRUE,
       xlab = "Prevalence of Lameness", smlab = "",
       squaresize = 0.5, col.square = "navy",
       col.diamond = "maroon", col.diamond.lines = "maroon",
       fs.hetstat = 10,
       print.Q = TRUE,
       print.pval.Q = TRUE,
       sortvar = seTE)
##
dev.off()