library(metafor) library(meta) library(tidyverse) prevalence_2020_cow_nomv = read.csv("prev_cow.cvs") prevalence_2020_cow_nomv ies.da = escalc(xi = nlameanimal, ni = ssizeanimal, data = prevalence_2020_cow_nomv, measure = "PFT", add = 0) pes.da = rma(yi, vi, data = ies.da, method = "DL") pes = predict(pes.da, transf = transf.ipft.hm, targ = list(ni = prevalence_2020_cow_nomv$ssizeanimal)) print(pes.da, digits = 3) confint(pes.da, digits = 3) print(pes, digits = 3) ies.da.lcmbi = ies.da %>% filter(is.na(lcmbi) == FALSE) subganal.lcmbi = rma(yi, vi, data = ies.da.lcmbi, mods = ~lcmbi, method = "DL") pes.da.lcm = rma(yi, vi, data = ies.da.lcmbi, subset = lcmbi=="Records", method = "DL") pes.da.records = rma(yi, vi, data = ies.da.lcmbi, subset = lcmbi=="Locomotion Scoring Method", method = "DL") pes.subg.lcmbi = predict(subganal.lcmbi, transf = transf.ipft.hm, targ = list(ni = prevalence_2020_cow_nomv$ssizeanimal)) dat.samevar = data.frame(estimate = c((pes.da.lcm$b)[1], (pes.da.records$b)[1]), stderror = c((pes.da.lcm$se)[1], (pes.da.records$se)[1]), tau2 = subganal.lcmbi$tau2) pes.da.lcmbi = rma(estimate, sei = stderror, method = "DL", data = dat.samevar) pes.lcmbi = predict(pes.da.lcmbi, transf = transf.ipft.hm, targ = list(ni = prevalence_2020_cow_nomv$ssizeanimal)) print(pes.da.lcm, digits = 3) # display subgroup 1 summary effect size print(pes.da.records, digits = 3) # display subgroup 2 summary effect size print(subganal.lcmbi, digits = 3) # display subgroup analysis results print(pes.lcmbi, digits = 3) # display recomputed summary effect size subganal.lcmbi = rma(yi, vi, data = ies.da.lcmbi, mods = ~lcmbi, method = "DL") pes.summary = metaprop(nlameanimal, ssizeanimal, author, data = ies.da.lcmbi, sm = "PFT", method.tau = "DL", method.ci = "NAsm", byvar = lcmbi, tau.common = TRUE, tau.preset = sqrt(subganal.lcmbi$tau2)) precision = sqrt(ies.da.lcmbi$vi) jpeg(file = "forest_cow_cow_lcmbi.jpeg", width = 11, height = 12, units = "in", res = 300) ## forest(pes.summary, xlim = c(0, 1), pscale = 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 = "", weight.study = "random", squaresize = 0.5, col.square = "navy", col.diamond = "maroon", col.diamond.lines = "maroon", pooled.totals = FALSE, comb.fixed = FALSE, fs.hetstat = 10, print.tau2 = TRUE, print.Q = TRUE, print.pval.Q = TRUE, print.I2 = TRUE, digits = 2, sortvar = precision) ## dev.off()