Dear group,

I am fairly new to Stan/Bayesian statistics and maybe my question is stupid. Aynways, here is my problem:

I use a partially-pooled model for binary repeated trials to estimate annual forest disturbance rates (i.e., proportion of forests that were disturbed by insects, wind, etc.). My data set has 32 rows (=32 years), one column indicating the ‘successes’ (=number of forested plots disturbed), and one column indicating the number of ‘trials’ (=total number of forested plots). The data is attached. For building the model I follow the suggestions by Bob Carpenter: http://mc-stan.org/users/documentation/case-studies/pool-binary-trials.html. The model ist also attached (see model_base.stan).

Up to here, everything is working nicely! However, I know that while interpreting the plots we sometimes omit a disturbances (because it was to ephemeral) and we sometimes comit a disturbances (labeling a plot as disturbed even though it is not). Hence, the number of ‘successes’ (=number of forested plots disturbed) is just an estimation of the unknown ‘real’ number. I would like to integrate this into my model to get a better representation of uncertainties.

My first try was to model the response (i.e., number of ‘successes’) itself as binomial distributed. Hence, instead of:

model {

mu ~ normal(-1, 1); // hyperprior

sigma ~ normal(0, 1); // hyperprior

alpha_std ~ normal(0, 1); // prior

y ~ binomial_logit(K, mu + sigma * alpha_std); // likelihood

}

I tried something like this:

model {

int y_obs[N];

mu ~ normal(-1, 1); // hyperprior

sigma ~ normal(0, 1); // hyperprior

alpha_std ~ normal(0, 1); // prior

for (i in 1:N)

y_obs[i] ~ binomial(K[i], y[i] / K[i]); // With y / K being the estimated success-rate

y_obs ~ binomial_logit(K, mu + sigma * alpha_std); // likelihood

}

This model, however, does not sample. I am stuck here as I a) do not understand why it does not work and b) how I should proceed. Any suggestions?

Here the code and data:

model_base.stan (885 Bytes)

forestplots.csv (441 Bytes)

disturbance_rate_model.R (800 Bytes)