# 2012.11.08

# data from Fig 2 in

# Emotional distress in infertile women and failure of assisted 
# reproductive technologies: meta-analysis of prospective 
# psychosocial studies

# Boivin J  BMJ  BMJ 2011;342:d223 doi:10.1136/bmj.d223


x=scan()
2.6 -0.40 -0.84  0.04
4.5	-0.17 -0.50  0.17
1.3	-0.34 -0.98  0.29
7.3	0.03 -0.24  0.29
1.0	-0.58 -1.31  0.15
2.3	0.04 -0.43  0.51
20.7	0.00 -0.16  0.16 
4.0	-0.23 -0.58  0.13
1.4	-0.43 -1.04  0.17
26.5	0.05 -0.09  0.19
18.7	-0.02 -0.19  0.15
2.4	0.32 -0.14  0.78
1.7	-0.18 -0.73  0.38
5.6	-0.23 -0.54  0.07

#
d=matrix(x,ncol=4,byrow=T)
colnames(d)=c("weight", "theta.hat50", "theta.hat2.5", "theta.hat97.5")
ds=data.frame(d)

ds$theta.hat.var = ( (ds$theta.hat97.5 - ds$theta.hat2.5)/(2*1.96) )^2 ;

ds$info = 1 / ds$theta.hat.var ; ds


pooled.theta.hat50 = sum(ds$theta.hat50*ds$info) / sum( ds$info )

info.pooled.theta.hat50 =  sum( ds$info )  

se.info.pooled.theta.hat50 =  sqrt(1/info.pooled.theta.hat50); 
se.info.pooled.theta.hat50

round(pooled.theta.hat50 + c(-1,0,1) *1.96 * se.info.pooled.theta.hat50 , 2)



# Funnel plot


se.max = max(sqrt(1/ds$info))

plot( ds$theta.hat50,  sqrt(1/ds$info), ylim=c(0, se.max),
  xlim=c(-0.6,0.6) , ylab="SE" )

segments( -1.96*se.max , se.max , 0,0 )
segments(  1.96*se.max , se.max , 0,0 )
segments(     0*se.max , se.max , 0,0 , col="grey70")





##### MCMC random effects model...


MODEL=readLines(,14)
model
{
  for( i in 1:14) {
       theta.hat.50[i] ~ dnorm(mu.theta[i],inv.var.error[i])
       inv.var.error[i]  <- 1/theta.hat.var[i]
       mu.theta[i] <- mu.overall + alpha[i] 
       alpha[i] ~ dnorm(0, inv.var.between.studies)
  }     
  mu.overall ~ dnorm(0, 0.001)

  sigma.between.studies  ~ dunif(0,1000)
  var.between.studies      <- sigma.between.studies  * sigma.between.studies 
  inv.var.between.studies  <- 1/var.between.studies
}

writeLines(MODEL,'MetaAnalysisModel.txt')

DATA=list(theta.hat.50 = ds$theta.hat50, theta.hat.var = ds$theta.hat.var )

library('rjags')

jags <- jags.model('MetaAnalysisModel.txt',
                   data = DATA,
                   inits=list(alpha = c(-5,-5,-5,-5,-5,-5,-5, 5,5,5,5,5,5,5), mu.overall= 1),
                   n.chains = 1,
                   n.adapt = 1000)
 
update(jags, 40000)
 
results=jags.samples(jags, 
  c('mu.overall', 'var.between.studies', 'sigma.between.studies', 'alpha'), 1000)
str(results)

summary(results$mu.overall[1,,1]) 
hist(results$mu.overall[1,,1],breaks=seq(-0.3,0.3,0.01)) 
 

round(mean(results$mu.overall[1,,1])+c(-1,0,1)*1.96*sqrt(var( results$mu.overall[1,,1])),2)

summary(results$var.between.studies[1,,1])
summary(results$sigma.between.studies[1,,1])
                      
          
# shrinkage  ( cf Efrom and Morris )

plot(c(-0.6, 0.6), c(0,1.5), col="white", 
  ylab="", xlab="Revised", yaxt="n",
frame=F)
text(0.6,0.98,"observed", adj=c(1,1))
text(-0.6,1.1,"SE", adj=c(0,0), srt=90)
segments(-0.6,1,0.6,1)
theta.hat.posterior=rep(NA,14)
COL=rainbow(14)
for (i in 1:14) {
 theta.hat.posterior[i]= mean(results$alpha[i,,1])
 segments(ds$theta.hat50[i], 1, theta.hat.posterior[i], 0 , col=COL[i])
 points(c(ds$theta.hat50[i]), c(1), pch=19 , cex=0.75, col=COL[i])
 se = sqrt( ds$theta.hat.var[i] )
 segments(ds$theta.hat50[i], 1, ds$theta.hat50[i], 1+ se, , col=COL[i] )
 }

#100.00	-0.04 -0.11  0.03