#
# Confidence Interval for random-effects meta analysis
# based on the Bartlett corrected likelihood ratio
# Noma H. (Statistics in Medicine 2011, 30: 3304-3312)  
#
noma<-function(y, var, mu2, t2, clow, cup){
V<-var
alpha<-0.95
a <- 1 - (1 - alpha) / 2
k <- length(y)
R <- qchisq(alpha,df=1)
ci2<-numeric(2)
ci2[1]<-clow
ci2[2]<-cup
#
mle<-function(y, V, mu, t2){
 t5 <- t2 + 0.01
 par0 <- c(mu, t5)
 for(h in 1:2000){
   a1 <- ( (y - mu) * (y - mu) - V ) / ( (V + t5) * (V + t5) )
   a2 <- 1 / ( (V + t5) * (V + t5) )
   b1 <- y / (V + t5)
   b2 <- 1 / (V + t5)
   mu <- sum(b1) / sum(b2)
   t5 <- sum(a1) / sum(a2)
   t5 <- max(0, t5)
   par1 <- c(mu, t5)
   dif1 <- max( abs(par1 - par0)  )
   if(dif1 < 10^-12) break
   par0 <- par1
  }
 t5 <- max(0, t5)
 return(c(mu, t5))
}
#
rmle<-function(y, V, mu, t2){
 t5 <- t2 + 0.01
 for(h in 1:2000){
   a1 <- ( (y - mu) * (y - mu) - V ) / ( (V + t5) * (V + t5) )
   a2 <- 1 / ( (V + t5) * (V + t5) )
   t6 <- sum(a1) / sum(a2)
   dif1 <- abs(t6 - t5) / abs(t6)
   if(dif1 < 10^-12) break
   t5 <- t6
  }
 t5 <- max(0, t5)
 return(t5)
}
#
# Likelihood ratio ( bisectional method )
#
ml <- mle(y, V, mu2, t2)   # Maximum likelihood estimate 
mlike0 <- sum( log( dnorm(y, mean = ml[1], sd = sqrt( V + ml[2] ) ) ) )
mu3 <- ml[1]
t3  <- ml[2]
ci3 <- numeric(2)

c0 <- ci2[1] - 5
c1 <- mu3
for(g in 1:400){
    c2 <- (c0 + c1) / 2
    t4 <- rmle(y,V,c2,t3)
    mlike1 <- sum( log( dnorm(y, mean = c2, sd = sqrt( V + t4 ) ) ) )
    LR1 <- -2*(mlike1 - mlike0)
    if( LR1 >= R ) c0 <- c2
    if( LR1 < R) c1 <- c2
    if( abs(LR1 - R) < 10^-12) break    
    if(g == 400) c2 <- NaN
  }
ci3[1] <- c2

c0 <- mu3
c1 <- ci2[2] + 5
for(g in 1:400){
    c2 <- (c0 + c1) / 2
    t4 <- rmle(y,V,c2,t3)
    mlike1 <- sum( log( dnorm(y, mean = c2, sd = sqrt( V + t4 ) ) ) )
    LR2 <- -2*(mlike1 - mlike0)
    if( LR2 >= R ) c1 <- c2
    if( LR2 < R) c0 <- c2
    if( abs(LR2 - R) < 10^-8) break    
    if(g == 400) c2 <- NaN
  }
ci3[2] <- c2
#
# Bartlett correction 
#
ci5 <- numeric(2)
c0 <- ci3[1] - 5
c1 <- mu3

for(g in 1:400){
    c2 <- (c0 + c1) / 2
    t4 <- rmle(y,V,c2,t3)
    mlike1 <- sum( log( dnorm(y, mean = c2, sd = sqrt( V + t4 ) ) ) )
    LR0 <- -2*(mlike1 - mlike0)
    u1 <- (V + t4)^-1
    u2 <- (V + t4)^-2
    u3 <- (V + t4)^-3
    U <- 1 + 2 * sum(u3) / ( sum(u1) * sum(u2) )
    LR1 <- LR0 / U 
    if( LR1 >= R ) c0 <- c2
    if( LR1 < R) c1 <- c2
    if( abs(LR1 - R) < 10^-12) break    
    if(g == 400) c2 <- NaN
  }
ci5[1] <- c2

c0 <- mu3
c1 <- ci3[2] + 5
for(g in 1:400){
    c2 <- (c0 + c1) / 2
    t4 <- rmle(y,V,c2,t3)
    mlike1 <- sum( log( dnorm(y, mean = c2, sd = sqrt( V + t4 ) ) ) )
    LR0 <- -2*(mlike1 - mlike0)
    u1 <- (V + t4)^-1
    u2 <- (V + t4)^-2
    u3 <- (V + t4)^-3
    U <- 1 + 2 * sum(u3) / ( sum(u1) * sum(u2) )
    LR2 <- LR0 / U 
    if( LR2 >= R ) c1 <- c2
    if( LR2 < R) c0 <- c2
    if( abs(LR2 - R) < 10^-8) break    
    if(g == 400) c2 <- NaN
  }
ci5[2] <- c2
return( c(mu3, ci5[1], ci5[2]) )
}