##################################################
#Lipsey, M.W. & Wilson, D.B. (2001). 
#Computational techniques for meta-analysis data. 
#In Practical meta-analysis (pp.129--145). 
#Thousand Oaks, CA: Sage Publications.
##################################################

#----------------------------------------------#
#Step1: The mean, confidence interval, 
#           and homogeneity test (pp.129--133)
#----------------------------------------------#

##(a) dataset
es <- c(-.33, .32, .39,.31,.17,.64,-.33,.15,-.02, 0)   #Effect Size
w  <- c(11.905,28.571,58.824,29.412,13.889,8.547,9.804,10.753,83.333,14.925)
k  <- length(es)                                       #number of study


##(b) Table 7.1
v    <- 1/w
wes  <- w * es
wessq <- w * es^2
wsq  <- w^2
random    <- factor(c(0,0,0,0,0,0,1,1,1,1))
intensity <- c(7,3,7,5,7,7,4,4,5,6)
data.frame(es, v, w, wes, wessq, wsq, random, intensity)


##(c) sum of the variables w, wes, wessq, and wsq.
sum.w     <- sum(w)
sum.wes   <- sum(wes)
sum.wessq <- sum(wessq)
sum.wsq   <- sum(wsq)


##(d) weighted mean effect size and standard error
pooled.es  <- sum.wes/sum.w       #weighted mean effect size
pooled.se  <- sqrt(1/sum.w)       #standard error of the mean effect size
z.value    <- abs(pooled.es)/pooled.se #z-test
lower.p.es <- pooled.es + qnorm(0.05/2, lower.tail=T) * pooled.se #95% confidence interval
upper.p.es <- pooled.es + qnorm(0.05/2, lower.tail=F) * pooled.se
p.value    <- pnorm(z.value, lower.tail=F) * 2


##(e) homogeneity test
Q         <- sum.wessq - sum.wes^2/sum.w
p.value.Q <- pchisq(Q, k-1, lower.tail=F)


##(f) 関数 library(rmeta),library(meta)を使用 
library(rmeta)
library(meta)
#----------------------------------------------------------#
#meta.summaries(d, se, method="fixed")  library(rmeta)
#metagen(TE, seTE)                      library(meta)
#d,TE には各研究の効果量. 
#se, seTEには各研究の効果量の標準誤差を指定する
#ここでは，w = 1/se^2 (3.24式) から，標準誤差を求めている 
#----------------------------------------------------------#
meta1 <- meta.summaries(d=es, se=sqrt(1/w), method="fixed") #library(rmeta)
meta1
summary(meta1)
metagen(TE=es, seTE=sqrt(1/w))                              #library(meta)




#----------------------------------------------#
#Step2: Random effect model (pp.134--135)
#----------------------------------------------#

##(a) variance component
var.comp <- (Q - (k-1)) / (sum.w-sum.wsq/sum.w) #variance component
v.asta   <- v + var.comp


##(b) sum of the variables w, wes, wessq
w2         <- 1/v.asta
sum.w2     <- sum(w2)
sum.wes2   <- sum(w2*es)
sum.wessq2 <- sum(w2*es^2)


##(c) weighted mean effect size and standard error
pooled.es2  <- sum.wes2/sum.w2       #weighted mean effect size
pooled.se2  <- sqrt(1/sum.w2)       #standard error of the mean effect size
z.value2    <- abs(pooled.es2)/pooled.se2 #z-test
lower.p.es2 <- pooled.es2 + qnorm(0.05/2, lower.tail=T) * pooled.se2 #95% confidence interval
upper.p.es2 <- pooled.es2 + qnorm(0.05/2, lower.tail=F) * pooled.se2
p.value2    <- pnorm(z.value2, lower.tail=F) * 2


##(d) 関数 library(rmeta),library(meta)を使用 
meta2 <- meta.summaries(d=es, se=sqrt(1/w), method="random")
meta2
summary(meta2)
metagen(TE=es, seTE=sqrt(1/w))



#----------------------------------------------#
#Step3: Analog to ANOVA (pp.135--138)
#----------------------------------------------#

##(a) homogeneity Q for each group
tapply(w*es^2, random, sum)     #sum(wessq)
tapply(w*es, random, sum)^2     #sum(wes)^2
tapply(w, random, sum)          #sum(w)

Qj  <- tapply(w*es^2, random, sum)  - tapply(w*es, random, sum)^2/tapply(w, random, sum)


##(b) homogeneity test
j   <- nlevels(random)          #number of categories

Qw  <- sum(Qj)
df.qw <- k - j
p.value.qw <- pchisq(Qw, df.qw, lower.tail=F)

Qb    <- Q - Qw
df.qb <- j - 1
p.value.qb <- pchisq(Qb, df.qb, lower.tail=F)


##(c) weighted mean effect size for each group
pooled.es.j  <- tapply(w*es, random, sum) / tapply(w, random, sum)
pooled.se.j  <- sqrt(1/tapply(w, random, sum))
z.value.j    <- abs(pooled.es.j)/pooled.se.j
lower.p.es.j <- pooled.es.j + qnorm(0.05/2, lower.tail=T) * pooled.se.j #95% confidence interval
upper.p.es.j <- pooled.es.j + qnorm(0.05/2, lower.tail=F) * pooled.se.j
p.value.j    <- pnorm(z.value.j, lower.tail=F) * 2


##(d) 関数 library(meta)を使用 
meta3 <- metagen(TE=es, seTE=sqrt(1/w))
summary.meta(meta3, byvar=as.numeric(random), bylab=levels(random))



#----------------------------------------------#
#Step4: Weighted regression analysis (pp.138--140)
#----------------------------------------------#

##(a) least square method
reg1 <- lm(es~random+intensity, weights=w)
anova(reg1)
summary(reg1)


##(b) statistical test
se.b <- coef(summary(reg1))[,2]
mse  <- anova(reg1)[3,3]
B    <- coef(summary(reg1))[,1]
se.b <- se.b/sqrt(mse)
z.value.b    <- B/se.b
p.value.b <- pnorm(abs(z.value.b), lower.tail=F) * 2



#----------------------------------------------#
#Step5: graphing techniques (pp.142--145)
#----------------------------------------------#

##(a) histogram, stem-and-leaf, scatterplot, boxplot
par(mfrow=c(2,2))
hist(es)
stem(es)
plot(w~es)
boxplot(es~random)


##(b) funnel plot
funnelplot(x=es, se=sqrt(1/w))
X11()


##(c) error-bar chart
par(mfrow=c(1,2))
metaplot(mn=es, se=sqrt(1/w))             #library(rmeta)

meta3 <- metagen(TE=es, seTE=sqrt(1/w))
plot(meta3, xlab="Effect Size", comb.f=T) #library(meta)

graphics.off()