Meta-analysis with use of different priors, a simulation

Modeling
1

Hi,
I am new to stan and I am trying to to replicate the jags code below, that is for random eff. bayesian meta-analysis, that tests 2 different priors for the same parameter (for precision).

my data is
y.dat=c(0.33, 0.38, 0.22, 0.22, 0.23, -0.12 , 0.11) //mean
V=c(0.039, 0.041, 0.020, 0.065, 0.035, 0.009, 0.002) //variance
std=sqrt(V)
N=7 // number of studies in meta-analysis

my jags code is
cat("data{
# nstud=number of studies
# nprior=number of different priors for variance
# create replicates of datasets
for(j in 1:N){
for(k in 1:3){
y[j,k] <- y.dat[j]
P[j,k] <- 1/V[j]
}
}
}
model
{
# loop over number of priors k
for(k in 1:3){
# loop over number of studies j
for(j in 1:N){
y[j,k] ~ dnorm(delta[j,k],P[j,k])
delta[j,k]~ dnorm(d[k],prec[k])
}
### Define the priors
# prior for pooled effect
d[k]~ dnorm(0,0.0001)
}
# priors for precision
# prior 1 - dunif( 0, 10 ) on sigma
sigma[1] ~ dunif( 0, 10 )
sigma.sq[1] <- sigma[1] * sigma[1]
prec[1] <- 1/sigma.sq[1]

prior 2 - dunif( 0, 100 ) on sigma

sigma[2] ~ dunif( 0, 100 )  
sigma.sq[2] <- sigma[2] * sigma[2]
prec[2] <- 1/sigma.sq[2]
}", file="aspirinRE.txt")

my stan code is provided below, but it generates results that are way different from those I am getting from jags…I would greatly appreciate any help explaining what I did wrong

yy <- cbind(y.dat, y.dat) // I am replicating database for each prior in R
sstd = cbind(std, std)

stan.dat <- list(N = 7, y = yy, std = sstd)

// saved as MA_mean.stan
data {
int<lower=0> N; // number of trials data$N
matrix[N,2] y.dat; // estimated y, creating a replicate database for each prior
matrix[N,2] std; // std, creating a replicate database for each prior
}
parameters {
real d[2];
real<lower=0> sigma[2];
matrix[N,2] delta;
}

model {
for(k in 1:2){// looping over priors
for(j in 1:N){// looping over the studies
y.dat[j,k] ~ normal(delta[j,k],std[j,k]);
delta[j,k]~ normal(d[k],sigma[k]);
}

  d[k]~ normal(0,0.0001); // prior for the pooled effect
}
 sigma[1] ~ uniform(0,100);  //priors for sigmas
 sigma[2] ~ uniform( 0,10);

}
fit <- stan(file = ‘test3.stan’, data = stan.dat, iter = 5000, chains =3)

2

the issue seems resolved, the prior for the pooled effect should be modeled as

d[k]~ normal(0,100);
that is corresponding to JAGS version of d[k]~ dnorm(0,0.0001);

1 Like