# Initialization between (-2, 2) failed after 100 attempts.

# Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”

[1] “error occurred during calling the sampler; sampling not done”

Winbugs original Code

model{ # *** PROGRAM STARTS

for(i in 1:ns){ # LOOP THROUGH STUDIES

delta[i,1] <- 0 # treatment effect is zero for control arm

mu[i] ~ dnorm(0,.0001) # vague priors for all trial baselines

for (k in 1:2) { # LOOP THROUGH ARMS

r[i,k] ~ dbin(p[i,k],n[i,k]) # binomial likelihood

logit(p[i,k]) <- mu[i] + delta[i,k] # model for linear predictor

rhat[i,k] <- p[i,k] * n[i,k] # expected value of the numerators

dev[i,k] <- 2 * (r[i,k] * (log(r[i,k])-log(rhat[i,k])) #Deviance contribution

- (n[i,k]-r[i,k]) * (log(n[i,k]-r[i,k]) - log(n[i,k]-rhat[i,k])))

}

resdev[i] <- sum(dev[i,]) # summed residual deviance contribution for this trial

delta[i,2] ~ dnorm(d[2],tau) # trial-specific LOR distributions

}

totresdev <- sum(resdev[]) #Total Residual Deviance

d[1]<- 0 # treatment effect is zero for reference treatment

d[2] ~ dnorm(0,.0001) # vague prior for treatment effect

sd ~ dunif(0,5) # vague prior for between-trial SD

tau <- pow(sd,-2) # between-trial precision = (1/between-trial variance)

}

R stan Code:-

rm(list=ls())

d1 <- c(3, 14, 11, 127, 27, 6, 152, 48, 37, 188, 52, 47, 16, 45, 31, 38, 12, 6, 3, 40, 43, 39)

n1 <- c(39, 116, 93, 1520, 365, 52, 939, 471, 282, 1921, 583, 266, 293, 883, 147, 213, 122, 154, 134, 218,

364, 674)

d2 <- c(3, 7, 5, 102, 28, 4, 98, 60, 25, 138, 64, 45, 9, 57, 25, 33, 28, 8, 6, 32, 27, 22)

n2 <- c(38, 114, 69, 1533, 355, 59, 945, 632, 278, 1916, 873, 263, 291, 858, 154, 207, 251, 151, 174, 209, 391, 680)

deaths <- matrix(data = c(d1,d2),byrow = F,nrow = 22, ncol = 2)

patients <-matrix(data = c(n1,n2),byrow = F,nrow = 22, ncol = 2)

library(“rstan”)

set.seed(110)

JJ <- nrow(patients)

stan.dat <- list(J = JJ, K1=2, r = deaths, n = patients)

bayesian_model<- "data {

int<lower=0> J; // number of trials

int<lower=0> K1; // number of treatments

int r[J, K1] ; // deaths of patients in trial j in arm k

int n[J, K1] ; // se of log relative risk

}

parameters {

vector[J] mu;

real death_dist;

real<lower=0, upper=5> sdd;

}

model {

matrix[J, K1] p;

matrix[J, K1] delta;

for (m in 1:J)

delta[m,1]<- 0;

for (a in 1:J) for (b in 1:K1)

p[a,b] <- inv_logit(mu[a] + delta[a,b]);

death_dist ~ normal(0,1);

sdd ~ uniform(0,5);

for (s in 1:J)

mu[s] ~ normal(0, 1);

for (t in 1:J) for (u in 1:K1)

r[t,u] ~ binomial(n[t,u], p[t,u]);

for (x in 1:J)

delta[x,2]~ normal(death_dist,pow(sdd,-2));

}"

fit <- stan(model_code = bayesian_model, data = stan.dat, iter = 2000, chains = 4)

post <- extract(fit, permuted = TRUE)

quantile(exp(post$mu), probs = c(.025, .5, .975))

Please help me resolve this problem.