#1.2A simple simulation on Multicenter Studies for Meta-analysis

# 1.2 Data Simulatio
# 1.2.1.2 Data Generation and Manipulation
### data seting
#Number of Participants each arm
n =100
# Mean blood pressure at baseline
mu = 100
# standard deviations for blood pressure
sd =20
# Mean changes for blood pressure
mu.d = 10
# Mean age for participants
age.mu = 50
# sd of age for participants
age.sd = 10

###simulation for CTRL
#Fix the seed for random number generation
set.seed(123)
# Use "rnorm" to generate random normal
#simuate age
age = rnorm(n,age.mu,age.sd)
# simulate `baseline` blood pressure
bp.base = rnorm(n,mu,sd)
# simulate `endpoint` blodd pressure
bp.end = rnorm(n,mu,sd)
# Take the difference between endpoint and baseline
bp.diff =bp.end -bp.base
# put the data together using "cbind" to column-bind
dat4CTRL = round(cbind(age,bp.base,bp.end,bp.diff))
# print the first 6 observations
head(dat4CTRL)

###simualtion for Drug
#simuate age
age = rnorm(n,age.mu,age.sd)
# simulate `baseline` blood pressure
bp.base = rnorm(n,mu,sd)
# simulate `endpoint` blodd pressure
bp.end = rnorm(n,mu-mu.d,sd)
# THe change in blood pressure
bp.diff =bp.end -bp.base
# Make the data matrix
dat4drug = round(cbind(age,bp.base,bp.end,bp.diff))
# print the first 6 observations
head(dat4drug)

# Make a dataframe to hold all data
dat1 = data.frame(rbind(dat4CTRL,dat4drug))
#Make "TRT" as a factor for treatment
dat1$TRT = as.factor(rep(c("CTRL","Drug"),each = n))
#Make a "Center" to represent the center number
dat1$Center = 1

###check the dataset
# check the data dimension
dim(dat1)

#print the first 6 observations to see the variable names
head(dat1)

##this process of data generation into a function
data.generator =function(n,age.mu,age.sd,mu,mu.d,sd,center){
  #Date from CTRL
  age = rnorm(n,age.mu,age.sd)
  bp.base = rnorm(n,mu,sd)
  bp.end = rnorm(n,mu,sd)
  bp.diff = bp.end -bp.base
  dat4CTRL = round(cbind(age,bp.base,bp.end,bp.diff))
  #Date from Drug
  age = rnorm(n,age.mu,age.sd)
  bp.base = rnorm(n,mu,sd)
  bp.end = rnorm(n,mu-mu.d,sd)
  bp.diff = bp.end -bp.base
  dat4drug = round(cbind(age,bp.base,bp.end,bp.diff))
  #Put both data matrice\tilde{}s together
  dat = data.frame(rbind(dat4CTRL,dat4drug))
  #Make "TRT" as a factor for treatment
  dat$TRT = as.factor(rep(c("CTRL","Drug"),each =n))
  #Make a "Center" to represent the center number
  dat$Center = center
  # Return the simulated data
  dat
  #end of function
}

##simulation data
#center1
d1 =data.generator(n,age.mu,age.sd,mu,mu.d,sd,1)
#center 2
mu.d2 =13
d2 =data.generator(n,age.mu,age.sd,mu,mu.d2,sd,2)
#center 3
mu.d3 =15
d3 =data.generator(n,age.mu,age.sd,mu,mu.d3,sd,3)
#center 4
mu.d4 =8
d4 =data.generator(n,age.mu,age.sd,mu,mu.d4,sd,4)
#center 5
mu.d5 =10
d5 =data.generator(n,age.mu,age.sd,mu,mu.d5,sd,5)

#putting these data from 5 centers together
dat = data.frame(rbind(d1,d2,d3,d4,d5))
#Change `Center` from numeric to factor
dat$Center = as.factor(dat$Center)

#1.2.1.3 Basic R Graphics
# the normally distribution check
#call boxplot
boxplot(dat4CTRL,las =1,main = "Control Drug")
boxplot(dat4drug,las =1,main = "New Drug")

#the relationship between the blood pressure difference as a function of age 
# for each treatment to assess whether there exists a statistically 
# significant relationship in addition to a treatment difference
library(lattice)
#call xyplot function and print it
print(xyplot(bp.diff~age|Center*TRT,data=dat,xlab ="Age",strip = strip.custom(bg="yellow"),ylab="Blood Pressure Difference",lwd=3,cex=1.3,pch=20,type=c("p","r")))

# illustrate the treatment effect by center
# call bwplot
print(bwplot(bp.diff~TRT|Center,data=dat,xlab ="TRT",strip = strip.custom(bg="yellow"),ylab="Blood Pressure Difference",lwd=3,cex=1.3,pch=20,type=c("p","r")))

# illustrate the treatment effect is to group the treatment by center
#call bwplot
print(bwplot(bp.diff~Center|TRT,data = dat,xlab = "Center",strip = strip.custom(bg="yellow"),ylab ="Blood Pressure Difference",lwd = 3,cex =1.3, pch =20, type =c("p","r")))

# Data Analysis from Each Center
# perform an analysis of variance
#Model for Center 1
m.c1 = aov(bp.diff~TRT,data = dat[dat$Center == 1,])
# Print the summary
summary(m.c1)
#Model for Center 2
m.c2 = aov(bp.diff~TRT,data = dat[dat$Center == 2,])
# Print the summary
summary(m.c2)
#Model for Center 3
m.c3 = aov(bp.diff~TRT,data = dat[dat$Center == 3,])
# Print the summary
summary(m.c3)
#Model for Center 4
m.c4 = aov(bp.diff~TRT,data = dat[dat$Center == 4,])
# Print the summary
summary(m.c4)
#Model for Center 5
m.c5 = aov(bp.diff~TRT,data = dat[dat$Center == 5,])
# Print the summary
summary(m.c5)

# Data Analysis with Pooled Data from Five Centers
# call `aov` to fit the 3-way model[anlaysis of variance(ANOVA)]
lm1= aov(bp.diff~TRT*Center*age,data =dat)
summary(lm1)
# call `aov` to fit the reduced model
lm2 = aov(bp.diff~TRT+Center,data = dat)
summary(lm2)

# 1.2.2.3 A Brief Introduction to Meta-Analysis
# Get the study sample size
ndat = aggregate(dat$bp.diff,list(Center=dat$Center,TRT=dat$TRT),length)

ndat

#calculate the means by study
mdat = aggregate(dat$bp.diff,list(Center=dat$Center,TRT=dat$TRT),mean)

mdat

#calculate the standard deviations
sddat = aggregate(dat$bp.diff,list(Center=dat$Center,TRT = dat$TRT),sd)
#Print the SDs
sddat

# calculate the effect-size(ES) using mean-difference
#Call the library
library(metafor)
# Calculate the ESs
esdat = escalc(measure="MD",
               n1i = ndat$x[ndat$TRT == "Drug"],
               n2i = ndat$x[ndat$TRT == "CTRL"],
               mli = mdat$x[mdat$TRT == "Drug"],
               m2i = mdat$x[mdat$TRT == "CTRL"],
               sdli = sddat$x[sddat$TRT == "Drug"],
               sd2i = sddat$x[sddat$TRT == "CTRL"],append= T)
rownames(esdat) =ndat$Study[ndat$TRT=="TRT"]
#Print the ES dataframe
esdat