Why am I getting negative means from a hierarchical clustering model fitted to binomial data using RStan

Question

Apologies in advance - I'm new to Bayesian statistical methods and may not get all the terminology correct.

I am trying to run some Bayesian analysis on survey data using the rstan package in R but I'm unsure of how to interpret my output. My aim is to generate a probability distribution for each level of data depending on the response, where the response has predictors variables and there are levels in the data.

The response data I'm inputting are binomial but the resulting probability distribution of the posterior has negative values.

My data are structured as follows:

• 3524 responses (N)
• 23 geographical areas sampled treated as levels (L) where each response is associated with an area (ll)
• each area is scored 1-10 depending on it's level of deprivation and this is treated as a predictor variable (x[,1])
• the age group of each responder is treated as a predictor variable (x[,2])
• there are two predictor variables (D)
• yes or no response to a question in the survey (y)

So for each responder I have a response (1 = yes, 0 = no), a geographical area (1-23), a level of deprivation (1-10) and an age group (1-8 where 1 is the youngest age group and 8 is the oldest).

I'm trying to fit a hierarchical logistical model as follows:

data {
  int<lower=1> D;                         
  int<lower=0> N;                         
  int<lower=1> L;                         
  array[N] int<lower=0, upper=1> y;       
  array[N] int<lower=1, upper=L> ll;      
  array[N] row_vector[D] x;               
}
parameters {
  array[D] real mu;                       
  array[D] real<lower=0> sigma;           
  array[L] vector[D] beta;                
}
model {
  for (d in 1:D) {
    mu[d] ~ normal(0, 100);
    for (l in 1:L) {
      beta[l, d] ~ normal(mu[d], sigma[d]);
    }
  }
  for (n in 1:N) {
    y[n] ~ bernoulli(inv_logit(x[n] * beta[ll[n]]));
  }
}

I then run the model using the following code:

fit = stan(model_code=bern.stan, data=list(y=y, N=N, x=x, D=D, L=L, ll=ll), iter=5000)

I can print the fit as follows:

print(fit, probs=c(0.1, 0.9))

The print shows that the mean for each beta is negative (e.g., beta[2,1] where 2 is the level and 1 is the predictor variable age is -0.17). When I extract these data to plot the distribution it is negative.

My question is how can the mean and credible intervals be negative when I'm dealing with binomial data and positive predictor variables?

Answer 1

Probabilities cannot be negative, but as said in the comment, this is a logistic regression model, and if you transformed them using inv_logit, as your model does, they would not go outside the bounds.