###########################################
#Kline, R. B. (2004).
#Replication and meta-analysis
#In Beyond significance testing: Reforming data analysis methods in behavioral research (pp.247-271).
#Washington, DC: American Psychological Association.
###########################################

#------------------------#
#Step1: (8.5)式の確認 p.260
#------------------------#
Q1 <- 10; Q2 <- 25; k <- 15; sig.level <- 0.05

qchisq(p= sig.level, df=(k-1), lower.tail=F)
pchisq(q= Q1, df=(k-1), lower.tail=F)
pchisq(q= Q2, df=(k-1), lower.tail=F)



#------------------------#
#Step2: Mini Meta-Analysis p.262-265
#------------------------#

##(a) Table 8.2 の基礎データを作成
tab8.2 <- matrix(c(1,2,3,4,5,6,7,8,9,
                     1,1,1,2,2,3,3,3,3,
                     25,15,40,20,35,30,15,20,45,
                     20,25,40,20,35,35,15,20,45,
                     -.05,.05,-.1,.55,.65,.9,.7,.85,.75),ncol=5)
tab8.2 <- data.frame(tab8.2)
names(tab8.2)  <- c("Study","Illness.chronicity", "Treatment", "Control", "g")
tab8.2$Illness.chronicity <- factor(tab8.2$Illness.chronicity, labels=c("Chronic","Intermediate","Acute"))
tab8.2


##(b) Table 4.5 (p.108) 効果量の標準誤差関数の定義
g.std <-function(g,n1,n2){
      dfw <- n1+n2- 2
      N   <- n1+n2
      sqrt(g^2/(2*dfw) + N/(n1*n2))
}


##(c) 8.2式の定義 (p.259)
attach(tab8.2)
   tab8.2$g.var <- (g.std(g=g,n1=Treatment,n2=Control))^2
   tab8.2$w     <- 1/tab8.2$g.var                        #8.2式
   tab8.2$wg    <- tab8.2$w * tab8.2$g
   tab8.2                                                #Table 8.2 完成
detach(tab8.2)


##(d) pooled effect size から等質性の検定まで
attach(tab8.2)
   (Mg <- sum(wg)/sum(w))                                  #8.1式
   (Smg <- sqrt(1/sum(w)))                                 #8.3式

   #8.4式
   (lower.Mg <- Mg + Smg * qnorm(sig.level/2, lower.tail=T))
   (upper.Mg <- Mg + Smg * qnorm(sig.level/2, lower.tail=F))
   
   (Q <- sum((g-Mg)^2/g.var))                               #8.5式
   
   k <- nrow(tab8.2)                                        #degree of freedom
   qchisq(p= sig.level, df=(k-1), lower.tail=F)             #critical value
   pchisq(q= Q, df=(k-1), lower.tail=F)                     #p.value
detach(tab8.2)


##(e) Table 8.3の計算
attach(tab8.2)

  #8.2式
  (Ww <- tapply(w, Illness.chronicity, sum))

  #8.1式
  (Mgw <- tapply(wg,Illness.chronicity, sum)/Ww)

  #8.3式
  (Smgw <- sqrt(1/Ww))
  
  #8.4式
   (lower.Mgw <- Mgw + Smgw * qnorm(sig.level/2, lower.tail=T))
   (upper.Mgw <- Mgw + Smgw * qnorm(sig.level/2, lower.tail=F))

   #8.5式
   Qw <- rep(NA, length(Mgw))
         for(i in 1:length(Mgw)){
               Qw[i] <- sum((g[Illness.chronicity==names(Mgw)[i]] - Mgw[i]) ^2 /
               g.var[Illness.chronicity==names(Mgw)[i]])
         }
   names(Qw) <- names(Mgw)
   (dfw <- tapply(w, Illness.chronicity, length)-1)
   pchisq(q= Qw, df=dfw, lower.tail=F)                     #p.value

   #between group homogeneity test statistics
   (Qb <- sum(Ww * (Mgw-Mg)^2))
   pchisq(q=Qb, df= nlevels(Illness.chronicity)-1, lower.tail=F)
   qchisq(p= sig.level, df=nlevels(Illness.chronicity)-1, lower.tail=F)
   
   #within group homogeneity test statistics
   Qw2 <- rep(NA, length(Mgw))
   for(i in 1:length(Mgw)){
      Qw2[i] <- sum(w[Illness.chronicity==names(Mgw)[i]]*
       (g[Illness.chronicity==names(Mgw)[i]]-Mgw[i])^2
       )
       }

   names(Qw2) <- names(Mgw)
   Qw2
   Qw                      #計算方法が異なるだけ
   qchisq(p=sig.level, df=dfw, lower.tail=F)    #critical value
   
   #Q = Qb + Qw
   sum(Qw) + Qb
   Q
detach(tab8.2)


#------------------------#
#Step3: library(rmeta) を使用
#------------------------#
library(rmeta)
attach(tab8.2)
   meta1 <- meta.summaries(d=g, se=sqrt(g.var), method="fixed")
   meta1
   summary(meta1)
   metaplot(mn=g, se=sqrt(g.var),labels=Study, xlab="effect size")

   graphics.off()

   meta.summaries(d=g, se=sqrt(g.var), method="fixed",subset= Illness.chronicity=="Chronic")
   meta.summaries(d=g, se=sqrt(g.var), method="fixed",subset= Illness.chronicity=="Intermediate")
   meta.summaries(d=g, se=sqrt(g.var), method="fixed",subset= Illness.chronicity=="Acute")
detach(tab8.2)



#------------------------#
#Step4: library(meta) を使用
#------------------------#
library(meta)
attach(tab8.2)
      metagen(TE=g, seTE=sqrt(g.var))
      metagen(TE=g, seTE=sqrt(g.var), subset= Illness.chronicity=="Chronic")
      metagen(TE=g, seTE=sqrt(g.var), subset= Illness.chronicity=="Intermediate")
      metagen(TE=g, seTE=sqrt(g.var), subset= Illness.chronicity=="Acute")
detach(tab8.2)