rstatisticsstanrstan

I kept getting error message: Stan model 'anon_model' does not contain samples


The components of my stan models are as follows:

// The input data is an array 'treatment', 'total' and 'dead' of length 'N'.
data {
  int<lower=1> N;
  int<lower=0, upper=1> treatment[N];
  int<lower=0> total_cases[N];
  int<lower=0> dead_cases[N];
}

// The parameters accepted by the model. Our model
// accepts two parameters 'alpha' and 'beta'.
parameters {
  real<lower=0> alpha;
  real<lower=0> beta;
}


// 
model {
  for (i in 1:N) {
   dead_cases[i] ~ binomial_logit(total_cases[i], alpha + beta * treatment[i]);
  }
  // dead_cases ~ binomial_logit(total_cases, alpha + beta * treatment);
}

//
generated quantities {
  int dead_cases_sim[N];
  for (i in 1:N){
    dead_cases_sim[i] = binomial_rng(total_cases[i], alpha + beta * treatment[i]);
  }
  // dead_cases_rep = binomial_rng(total_cases, alpha + beta * treatment);
}
num_of_cases <- list(
  N = nrow(num_of_cases),
  total_cases = as.vector(num_of_cases$total_cases),
  dead_cases = as.array(num_of_cases$dead_cases),
  treatment = as.array(num_of_cases$Treatment))
fit <- stan(
  file = 'stan_EHT.stan',
  data = num_of_cases,
  chains = 4,
  warmup = 2000,
  iter = 4000,
  cores = 7
)

Sorry for the long text and ta. Basically what I'm trying to do is compare the effectiveness of a treatment (using binary code 1) with its corresponding control group (using binary code 0). The statistic I use is the number of dead cases (out of the total cases), so I build the binomial models (with input n = total cases, p = alpha + beta * treatment), hoping that I can see the difference between treatment and control.

What I've experimented and concluded up to now is that:


Solution

  • I solved this question by asking my supervisor, in short, the problem lies in the parameterizations of alpha + beta * treatment[i], I used the logit transformation previously, but I didn't transform it back somehow, therefore the parameter p = alpha + beta * treatment[i] is rejected everytime (since it's <0 sometimes).