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.