# Chapter: meta-analysis

# http://www.robin-beaumont.co.uk/virtualclassroom/book2data/cummingp245.rdata 

# load file from url:
load(url("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/cummingp245.rdata"))
# alternatively
Save the file locally taking a note of where you put it then use the following code to load the file load(file=file.choose()) 

# cumming in book2data cummingp245
attach(cummingp245)
names(cummingp245)
# [1] "author" "year"   "mean"   "sd"     "number" "method"
levels(cummingp245$method)
# [1] "alloc"  "na"     "random"

# install.packages("metafor")
library(metafor)

res <- rma(measure ="MN", mi=mean, sdi=sd,  ni=number, method="DL", data=cummingp245)
summary(res)
forest(res, showweight=TRUE)
forest(res, showweight=TRUE, slab=paste(author, year, sep=",   "))

funnel(res)
labbe(res)  # does not work for outcome measure = MN?

###### 
# knha=logical specifying whether the method by Knapp and Hartung (2003) should be 
# used for adjusting test statistics and confidence intervals (default is FALSE).
# tau2 = optional numerical value to specify the amount of (residual) 
# heterogeneity in a random- or mixed-effects model (instead of estimating it). 
# Useful for sensitivity analyses (e.g., for plotting results as a function of tē). 
# When unspecified, the value of tē is estimated from the data.
res <- rma(measure ="MN", mi=mean, sdi=sd,  ni=number, method="DL", data=cummingp245, knha=TRUE)
summary(res)
confint(res)
res2 <- rma(measure ="MN", mi=mean, sdi=sd,  ni=number, method="DL", data=cummingp245, knha=FALSE)
summary(res2)
confint(res2)

## Accounting for heterogeneity
#### now subgroup analysis (factor moderator)
res_alloc <- rma(measure ="MN", mi=mean, sdi=sd,  ni=number, method="DL", data=cummingp245, 
subset=(method=="alloc"), slab=paste(author, year, sep=",   "))
# need to specify slab content in above as subset will not work below
summary(res_alloc)
forest(res_alloc, showweight=TRUE, main="studies using NONrandom allocation")
res_random <- rma(measure ="MN", mi=mean, sdi=sd,  ni=number, method="DL", data=cummingp245, 
subset=(method=="random"),  slab=paste(author, year, sep=",   "))
summary(res_random)
forest(res_random, showweight=TRUE, main="studies using random allocation")

# to compare the two
data_reduced <- cummingp245[cummingp245$method !="na",]
data_reduced <- droplevels(data_reduced) # necessary to remove the na level even with no cases in it

res2 <- rma(measure ="MN", mi=mean, sdi=sd,  ni=number, 
      method="DL", mods = ~method, data=data_reduced )
summary(res2)
forest(res2, order = order(data_reduced$method))

res3 <- rma(measure ="MN", mi=mean, sdi=sd,  ni=number, 
      method="DL", mods = ~method + year, data=data_reduced )
summary(res3)

######### year insignificant so go back to res2
plot(res2)
baujat(res2)

#Random effects options:
# method="DL" = DerSimonian-Laird estimator
# method="HE" = Hedges estimator
# method="HS" = Hunter-Schmidt estimator
# method="SJ" = Sidik-Jonkman estimator
# method="ML" = maximum-likelihood estimator
# method="REML" = restricted maximum-likelihood estimator
# method="EB" = empirical Bayes (= Paule-Mandel) estimator
# method="GENQ" = generalized Q-statistic estimator