The components of my stan models are as follows:
stan_EHT.stan
)// 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:
in the generated quantities
, if I include any of the parameters (namely alpha
and beta
), I will get the error message "Stan model 'anon_model' does not contain samples.";
if I remove the whole generated quantities
, the code also works, meaning that other part of the code is correct.
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).