#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#
# Introduction to meta-analysis #
# ... 05th December 2025        #
#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#~#

# When you join the Zoom meeting:
# - Open the same pdf copy of the
#   notes (starting on page 84)
# - Open Rstudio (latest version)
# - access the following web page:
#   http://www.casc-platforms.com/MetaR
# - Make sure you have the two 
#   data sets saved in an appropriate 
#   place on your computer...

# remind me to start recording at the
# beginning! We'll begin at 9.30am.

# To run a line of code:
# Cntrl + Enter
# Command + Enter

# Getting started - page 87 ####

# 1. setting the working directory
# can go through drop-down menu.
# Session > Set working directory > Choose
setwd("C:/Users/deanl/OneDrive - University College London/Intro to meta-analysis")

# 2. create and save script file.

# 3. read in the data.
list.files()
fatigue <- read.csv("datafatigue.csv")
OME <- read.csv("dataOME.csv")
# or 
OME = read.csv("dataOME.csv")

# 4. install relevant packages.
# metafor
# how to install?
install.packages("metafor")
# load package after installation
library(metafor)

# test to see if installed successfully, run:
escalc()
# do you get the same error, or does it say
# "function not found"

# Exercise 1 - page 90 ####

# everyone done :)

# View, access and reformat - page 90 ####

View(OME)
# capital V
names(OME)
head(OME)

# access / use certain rows / columns
fatigue[,1]
# within fatigue data, give me 1st column
fatigue[1,]
# within fatigue data, give me 1st row
# [] square brackets - access part of something (object)
fatigue[,"studyname"]
fatigue$studyname

fatigue[c(1,2,3),]
# give me first 3 rows in fatigue data

# last thing ... subsetting the data
# only non-anti-TNF studies ...
# two ways (notice two == signs) ...
fatigue.TNF <- fatigue[fatigue$anti.TNF == 1,]
# or ...
fatigue.TNF <- subset(fatigue, anti.TNF==1)
# DO THIS ABOVE...

# please note - there is a popular package
# for data manipulation ... tidyverse

# numeric outcome meta-analysis - page 93 ####

# what does meta-analysis involve? steps?
# calculate effect estimates and standard errors
fatigue.TNF.SMD <- escalc(measure = "SMD", 
  m1i = mean.bio,
  m2i = mean.control,
  sd1i = sd.bio,
  sd2i = sd.control,
  n1i = total.bio,
  n2i = total.control,
  data = fatigue.TNF)

# version in the notes, but better to do as above
fatigue.TNF.SMD <- escalc(measure = "SMD", 
  m1i = fatigue.TNF$mean.bio,
  m2i = fatigue.TNF$mean.control,
  sd1i = fatigue.TNF$sd.bio,
  sd2i = fatigue.TNF$sd.control,
  n1i = fatigue.TNF$total.bio,
  n2i = fatigue.TNF$total.control)

# these are the new columns just created...
fatigue.TNF.SMD$yi
fatigue.TNF.SMD$vi

# access help page
?escalc

# next step, run meta-analysis
rma(fatigue.TNF.SMD)
# what's happened? methods? conclusions?

# break until 11.30am - see you soon!

# save results
fatigue.TNF.res <- rma(fatigue.TNF.SMD)
round(fatigue.TNF.res$beta,2)

# alternatives...
# for fixed effect analysis
rma(fatigue.TNF.SMD, method = "FE")
# for other tau2 estimator
rma(fatigue.TNF.SMD, method = "DL")

# final step - forest plot.
forest(fatigue.TNF.res)
# and with author names...
forest(fatigue.TNF.res, 
    slab = fatigue.TNF.SMD$studyname)
dev.copy(device = png, "myforestplot.png",
      width = 500, height = 750)
dev.off()

# Exercise 2 - page 105 ####

# Solutions (provisionally) at 12.05pm

# part (a) hint. change this line from earlier
fatigue.TNF <- fatigue[fatigue$anti.TNF == 0,]

# Harith's solution:

fatigue.nonTNF <- subset(fatigue, anti.TNF==0)

fatigue.nonTNF.SMD <- escalc(measure = "SMD", 
   m1i = mean.bio,
   m2i = mean.control,
   sd1i = sd.bio,
   sd2i = sd.control,
   n1i = total.bio,
   n2i = total.control,
   data = fatigue.nonTNF)

fatigue.nonTNF.res <- rma(fatigue.nonTNF.SMD)
fatigue.nonTNF.res

forest(fatigue.nonTNF.res)

# binary outcome meta-analysis - page 99 ####

# effect sizes and standard errors
OME.RR <- escalc(measure = "RR",
  ai = antibiotics.event, # number of outcomes in group 1
  ci = control.event, # outcomes in group 2
  n1i = antibiotics.total, # sample size in group 1
  n2i = control.total, # sample size in group 2
  data = OME)

OME.RR$yi # RR? No... ln(RR)
exp(OME.RR$yi) # RR
OME.RR$vi # variance
sqrt(OME.RR$vi) # standard error

# do meta-analysis
OME.res <- rma(OME.RR)
OME.res

# pooled result
OME.res$b # ln(RR)
exp(OME.res$b) # RR
exp(OME.res$ci.lb) 
exp(OME.res$ci.ub) 

# last step ... forest plot
forest(OME.res, 
    slab = paste(OME$author, "(", OME$year,")"),
    atransf = exp,
    at = log(c(0.05,.25,1,4,20)))


# exploring heterogeneity - page 106 ####

# sub group analysis - only categorical
# meta-regression - numeric / categorical

# meta-regression ####

fatigue.TNF.SMD$age

metareg <- 
  rma(fatigue.TNF.SMD, mods = ~ fatigue.TNF.SMD$age)
metareg
# equation...
# SMD = -0.3157 - 0.0020 * Age
# SMD = -0.3157 - 0.0020 * 0 = -0.3157
# SMD = -0.3157 - 0.0020 * 50 = -0.4157

# bubble plot
wi <- 1/fatigue.TNF.SMD$vi
size <- (wi+5)/25
plot(x = fatigue.TNF.SMD$age,
     y = fatigue.TNF.SMD$yi,
     cex = size)
abline(a = -0.3157, b = -0.0020)

# sub group analysis ####

# grouping the studies together.
OME$duration.3mth

# step 1 - meta-analysis of each sub group
subgroup0 <- rma(OME.RR[OME.RR$duration.3mth == 0,])
subgroup0
subgroup1 <- rma(OME.RR[OME.RR$duration.3mth == 1,])
subgroup1
# step 2 - test for differences between them
# two effect sizes
EffSize <- c(subgroup0$b, subgroup1$b)
EffErrors <- c(subgroup0$se, subgroup1$se)
rma(yi = EffSize, sei = EffErrors)

# is rma like any type of model?
AIC(metareg)