I am trying to replicating a random effect model from WinBUGS to RSTAN but getting error:-

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.

Stan’s error messages are great because they usually tell you how to resolve the problem. As it says “Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.”

So the random initialisation of your parameters that Stan uses by default isn’t working given the model you’ve specified and your parameter constraints. Try giving it some sensible starting values as a list.

I put some comments inline. There are some syntax-y problems going on here. Probly with the difference in expressing a model in WinBUGS and Stan.

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;

  // At this point delta[m, 2:K2] is full of NaNs
  // (unitialized variables in Stan are NaNs, variables defined in the model block are not parameters)
  for (a in 1:J)
    for (b in 1:K1)
      // This next line will just produce more NaNs
      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)
      // Because p is NaNs, this will make more NaNs
      r[t,u] ~ binomial(n[t,u], p[t,u]);

  for (x in 1:J)
    // normal in Stan is parameterized in terms of standard deviation
    delta[x, 2] ~ normal(death_dist, pow(sdd, -2));
}