library(metafor)
SRMA <- read.table(file="~/Dropbox/Lyme SRMA/R data/Final_Challenge_trials.txt",
                     sep = '\t', 
                     header = T, 
                     check.names = F)

#for lameness outcome
dat<- escalc(measure="OR", ai=lameness_cases_vax, ci=lameness_cases_ctrl, 
             n1i=total_vax,
              n2i=total_ctrl, data=SRMA, 
              #add=0.001, to="only0", 
             #drop00=T,
              append=TRUE)

#to calculate variance-covariance matrix accounting for correlated ORs 
dat$pti <-with(dat, lameness_cases_ctrl/total_ctrl)
dat$pci <-with(dat, lameness_cases_vax/total_vax)
dat$nli <-with(dat, total_vax)

calc.v <- function(x) {
  v <- matrix(1/(x$n1i[1]*x$pci[1]*(1-x$pci[1])), nrow=nrow(x), ncol=nrow(x))
  diag(v) <- x$vi
  v
}

V <- bldiag(lapply(split(dat, dat$study), calc.v))


# Check for NA's in the V matrix:
anyNA(V)
# If NA's exist, replace these with zero's for further calculations: 
V[is.na(V)] <- 0

res <-rma.mv(yi, V, random= list(~1|Publication_cluster, ~1|Within_study_cluster), data = dat)

#forest plot 
forest(res, slab = paste(dat$Authors),showweights=TRUE, atransf = exp,
       xlim = c(-16,8), at = log(c(0.01, 0.25, 1, 40)),
       ilab = cbind(dat$lameness_cases_vax, dat$total_vax, dat$lameness_cases_ctrl, 
                    dat$total_ctrl),
       ilab.xpos=c(-10.5,-9,-7.5,-6), cex=.75)

op <- par(cex= 0.75, font = 2)
text(-16, 21, "Author(s) and Year", pos = 4)
text(8, 21, "OR [95% CI]", pos = 2)
text(c(-10.5,-9,-7.5,-6), 21, c("Lame+", "Total", "Lame+", "Total"))
text(c(-9.75, -6.75), 22.5, c("Vaccinated", "Control"))
text(5.25, 21, "Weight")
text(-13.75, -2, "Overall I-squared= 37.1%, p=0.905")
par(op)


#calculating I2
W <- V
X <- model.matrix(res)
P <- W - W %*% X %*% solve(t(X) %*% W %*% X) %*% t(X) %*% W
100 * sum(res$sigma2) / (sum(res$sigma2) + (res$k-res$p)/sum(diag(P)))

funnel(res)

weights(res, type="diagonal")

#for subgroup analyses, but not needed because test for heterogeneity not significant
res<-rma.mv(yi, vi, data=dat, subset= (Commercially_available=="yes"),
            random=~1|Cluster_no)