#############################################
#Cohn L.D., & Becker, B.J. (2003).
#How meta-analysis increases statistical power.
#Psychological Methods, 8, 243--253.
#############################################

#------------------------------------------#
# Standardized Mean Differences p.245-246
#------------------------------------------#

##関数の定義
d.meta <- function(n1,n2, d, sig.level=.05,theta1=.36,theta=0,digits=2){
  v  <- ((n1+n2)/(n1*n2)) + (d^2/(2*(n1+n2)))
  w  <- 1/v
  
  T. <- cumsum(w*d) / cumsum(w)

  V. <- 1/cumsum(w)

  lower.T <- T. + sqrt(V.) * qnorm(sig.level/2, lower.tail=T)
  upper.T <- T. + sqrt(V.) * qnorm(sig.level/2, lower.tail=F)

  lamda <- (theta1-theta)/sqrt(V.)
  Power <- 1-pnorm(qnorm(sig.level/2,lower.tail=F)-lamda) + pnorm(qnorm(sig.level/2,lower.tail=T)-lamda)

  ##output
  result <- round(data.frame(d, n1,n2, N= cumsum(c(n1+n2)),ES=T., lower=lower.T, upper=upper.T, Power), digits)
  return(result)
}


##Table 1
n1 <- c(28,71,678,704,108,139,55,303,88,559,156,906)
n2 <- c(40,118,673,672,125,147,101,376,178,373,147,880)
d  <- c(.03,0.53,.67,.06,.33,-.08,.75,.48,.61,.22,.23,.43)
d.meta(n1=n1,n2=n2,d=d,theta1=.36,theta=0)


##Figure 1 A
fig1a <- d.meta(n1=n1,n2=n2,d=d,theta1=.36,theta=0)
attach(fig1a)
     plot(ES, ylim=c(-.6,.8),xlab="",xaxt="n",pch=18)
     arrows(x0=1:length(ES),y0=lower, x1=1:length(ES),y1=upper,code=3,angle=90, length=.02)
detach(fig1a)
graphics.off()



#------------------------------------------#
# Pearson Correlations p.246-248
#------------------------------------------#

##関数の定義
r.meta <- function(n, r, sig.level=.05,theta1=.25,theta=0,digits=2){
  z  <- atanh(r)
  w  <- n-3

  Zr. <- cumsum(w*z) / cumsum(w)
  V. <- 1/cumsum(w)

  lower.Zr. <- Zr. + sqrt(V.) * qnorm(sig.level/2, lower.tail=T)
  upper.Zr. <- Zr. + sqrt(V.) * qnorm(sig.level/2, lower.tail=F)

  lamda <- (theta1-theta)/sqrt(V.)
  Power <- 1-pnorm(qnorm(sig.level/2,lower.tail=F)-lamda) + pnorm(qnorm(sig.level/2,lower.tail=T)-lamda)

  ##output
  result <- round(data.frame(r, n, N= cumsum(n), ES=Zr., lower=lower.Zr., upper=upper.Zr., Power), digits)
  return(result)
}


##Table 2
n <- c(65,30,93,36,57,30,40,13,44,58,125,10,13)
r <- c(.17,-.45,-.47,-.13,-.24,-.15,-.25,-.66,-.25,-.23,-.18,.18,-.74)
r.meta(n=n,r=r,theta1=.25,theta=0)


##Figure 1 B
fig1b <- r.meta(n=n,r=r,theta1=.25,theta=0)
attach(fig1b)
     plot(ES, ylim=c(-.6,.6),xlab="",xaxt="n",pch=18)
     arrows(x0=1:length(ES),y0=lower, x1=1:length(ES),y1=upper,code=3,angle=90, length=.02)
detach(fig1b)
graphics.off()



#------------------------------------------#
# Odds Ratios p.248-249
#------------------------------------------#

##関数の定義
odds.meta <- function(n1,n2,n3,n4, odds, sig.level=.05,theta1=.11,theta=0,digits=3){
  v  <- 1/n1 + 1/n2 + 1/n3 + 1/n4
  w  <- 1/v

  L. <- cumsum(w*odds) / cumsum(w)
  V. <- 1/cumsum(w)

  lower.L. <- L. + sqrt(V.) * qnorm(sig.level/2, lower.tail=T)
  upper.L. <- L. + sqrt(V.) * qnorm(sig.level/2, lower.tail=F)

  lamda <- (theta1-theta)/sqrt(V.)
  Power <- 1-pnorm(qnorm(sig.level/2,lower.tail=F)-lamda) + pnorm(qnorm(sig.level/2,lower.tail=T)-lamda)

  ##output
  result <- round(data.frame(odds, n= n1+n2+n3+n4, N= cumsum(n1+n2+n3+n4), ES=L., lower=lower.L., upper=upper.L., Power), digits)
  return(result)
}


##Table 3
odds  <- c(.329, .385, .220, .222, .226, -.125, .111)
n1    <- c(566,714,730,285,725,2021,7017)
n2    <- c(49,44,102,32,85,246,1570)
n3    <- c(557,707,724,271,354,2038,6880)
n4    <- c(67,64,126,38,52,219,1720)

odds.meta(n1=n1,n2=n2,n3=n3,n4=n4, odds=odds,theta1=.11,theta=0)


##Figure 1 C
fig1c <- odds.meta(n1=n1,n2=n2,n3=n3,n4=n4, odds=odds,theta1=.11,theta=0)
attach(fig1c)
     plot(ES, ylim=c(-.1,.8),xlab="",xaxt="n",pch=18)
     arrows(x0=1:length(ES),y0=lower, x1=1:length(ES),y1=upper,code=3,angle=90, length=.02)
detach(fig1c)
graphics.off()



#------------------------------------------#
# Random-Effects Model p.249-251
#------------------------------------------#

##関数の定義
d.meta.random <- function(n1,n2, d, sig.level=.05,theta1=.36,theta=0,digits=2){

  v  <- ((n1+n2)/(n1*n2)) + (d^2/(2*(n1+n2)))
  w  <- 1/v

  T. <- cumsum(w*d) / cumsum(w)
  V. <- 1/cumsum(w)

  Q <-  cumsum(w*d^2)-(cumsum(w*d)^2)/cumsum(w)
  c  <- cumsum(w) - (cumsum(w^2)/cumsum(w))
  k  <- 1:length(d)
  tau.sq <- (Q - (k-1))/c
  
  lower.T <- upper.T <- lamda <- Power <- T.
  for(i in 2:length(d)){
    w  <- 1/(v[1:i] + tau.sq[i])
    T.[i] <- sum(w*d[1:i])/sum(w)
    V.[i] <- 1/sum(w)
    lower.T[i] <- T.[i] + sqrt(V.[i]) * qnorm(sig.level/2, lower.tail=T)
    upper.T[i] <- T.[i] + sqrt(V.[i]) * qnorm(sig.level/2, lower.tail=F)
    lamda[i] <- (theta1-theta)/sqrt(V.[i])
    Power[i] <- 1-pnorm(qnorm(sig.level/2,lower.tail=F)-lamda[i]) + pnorm(qnorm(sig.level/2,lower.tail=T)-lamda[i])    
  } 
  tau.sq[1] <- T.[1] <- lower.T[1] <- upper.T[1] <- Power[1] <- NA

  ##output
  result <- round(data.frame(d, n1,n2, N= cumsum(c(n1+n2)), var.comp =tau.sq ,ES=T., lower=lower.T, upper=upper.T, Power), digits)
  return(result)
}


##Table 1
n1 <- c(28,71,678,704,108,139,55,303,88,559,156,906)
n2 <- c(40,118,673,672,125,147,101,376,178,373,147,880)
d  <- c(.03,0.53,.67,.06,.33,-.08,.75,.48,.61,.22,.23,.43)
d.meta.random(n1=n1,n2=n2,d=d,theta1=.36,theta=0)