setwd( "C:/Users/kenrice/OneDrive - UW/Desktop/SISGBayes/exercises" )

# session 7 exercise, part 1

study.names <- c("AGES","ARIC","CHS","FHS","RS","RES")
betahat <- c( -0.37, -0.23, -0.31, -0.79, -0.13, -0.96) 
stderr  <- c(  0.26,  0.16,  0.33,  0.16,  0.25,  0.39)
n       <- c(  3219,  8047,  3277,  8096,  4737,  1760)

#install.packages("metafor")
library("metafor")
rma1 <- rma.uni(yi=betahat, vi=stderr^2, method="FE")
rma2 <- rma.uni(yi=betahat, vi=stderr^2, method="DL")

forest(rma1, slab=study.names, ylim=c(-2.5,8.5), xlim=c(-3,2))
addpoly(rma2, mlab="DSL")

# Bayes meta in Stan

cat(file="metaanalysis.stan", "
data {
  int p; //the number of studies
  vector[p] betahat; // study-specific estimates
  vector[p] sigma; // study-specific standard errors
  real T; // upper bound of uniform prior on tau
}
parameters {
  vector[p] beta; // study-specific parameters  
  real mu; // mean of the mixing distribution for beta[i]'s
  real tau; // SD of the mixing distribution
}
transformed parameters {
  vector[p] invsigma2;
  vector[p] betafpart;
  real betaf;
  for(i in 1:p){ invsigma2[i] = pow(sigma[i],-2);}
  for(i in 1:p){ betafpart[i] = beta[i]*invsigma2[i];}
  betaf = sum(betafpart[1:p])/sum(invsigma2[1:p]); // precision-weighted mean of study-specific effects
}
model {  
  for(i in 1:p){
	betahat[i] ~ normal(beta[i], sigma[i]);
	beta[i] ~ normal(mu, tau);
   }
   mu ~ normal(0,1);
   tau ~ uniform(0,T);
}
")

# do the MCMC, store the results
library("rstan")
stan.meta1 <- stan(file = "metaanalysis.stan", 
data = list(betahat=betahat, sigma=stderr, p=6, T=1),
iter = 100000, chains = 1, pars=c("betaf", "mu"))

print(stan.meta1)

cat(file="metaanalysisJAGS.txt", "
model {  
  for(i in 1:p){
	prec[i] <- 1/sigma[i]/sigma[i]
	betahat[i] ~ dnorm(beta[i], prec[i]);
	beta[i] ~ dnorm(mu, invtau);
	betafpart[i] <- beta[i]*prec[i]
   }
   betaf <- sum(betafpart[1:p])/sum(prec[1:p])
   invtau <- 1/tau
   mu ~ dnorm(0,1);
   tau ~ dunif(0,T);
}
")

library("rjags")
jags1 <- jags.model("metaanalysisJAGS.txt", 
	data = list(betahat=betahat, sigma=stderr, p=6, T=1/5) )
update(jags1, 10000)
summary( coda.samples(jags1, c("betaf","mu"), n.iter=10000) )
plot( coda.samples(jags1, c("betaf","mu"), n.iter=10000) )

post1 <- as.data.frame(coda.samples(jags1, c("betaf","mu"), n.iter=10000)[[1]])

forest(rma1, slab=study.names, ylim=c(-2.5,8.5), xlim=c(-3,2))
addpoly(x = mean( post1$betaf ) , sei = sd( post1$betaf ), mlab="Bayes FE", 
rows=-2)

forest(rma1, slab=study.names, ylim=c(-2.5,8.5), xlim=c(-3,2))
addpoly(x = mean( post1$mu ) , sei = sd( post1$mu ), mlab="Bayes RE", 
        rows=-2)