How can I model a proportion with BUGS/JAGS/STAN?

Question

I am trying to build a model where the response is a proportion (it is actually the share of votes a party gets in constituencies). Its distribution is not normal, so I decided to model it with a beta distribution. I also have several predictors.

However, I don't know how to write it in BUGS/JAGS/STAN (JAGS would be my best choice, but it doesn't really matter). My problem is that I make a sum of parameters by predictors, but then what can I do with it?

The code would be something like this (in JAGS syntax), but I don' know how to "link" the y_hat and y parameters.

for (i in 1:n) {
 y[i] ~ dbeta(alpha, beta)

 y_hat[i] <- a + b * x[i]
}

(y_hat is just the cross-product of parameters and predictors, hence the deterministic relationship. a and b are the coefficients which I try to estimate, x being a predictor).

Thanks for your suggestions!

Answer 1

The beta regression approach is to reparameterize in terms of $\mu$ and $\phi$. Where $\mu$ will be the equivalent to y_hat that you predict. In this parameterization you will have $\alpha=\mu\times\phi$ and $\beta=(1-\mu) \times \phi$. Then you can model $\mu$ as the logit of the linear combination. $\phi$ can either have its own prior (must be greater than 0), or can be modeled on covariates as well (choose a link function to keep it greater than 0, such as exponential).

Possibly something like:

for(i in 1:n) {
  y[i] ~ dbeta(alpha[i], beta[i])
  alpha[i] <- mu[i] * phi
  beta[i]  <- (1-mu[i]) * phi
  logit(mu[i]) <- a + b*x[i]
}
phi ~ dgamma(.1,.1)
a ~ dnorm(0,.001)
b ~ dnorm(0,.001)

Answer 2

Greg Snow gave a great answer. For completeness, here is the equivalent in Stan syntax. Although Stan has a beta distribution that you could use, it is faster to work out the logarithm of the beta density yourself because the constants log(y) and log(1-y) can be calculated once at the outset (rather than every time that y ~ beta(alpha,beta) would be called). By incrementing the reserved lp__ variable (see below), you can sum the logarithm of the beta density over the observations in your sample. I use the label "gamma" for the parameter vector in the linear predictor.

data {
  int<lower=1> N;
  int<lower=1> K;
  real<lower=0,upper=1> y[N];
  matrix[N,K] X;
}
transformed data {
  real log_y[N];
  real log_1my[N];
  for (i in 1:N) {
    log_y[i] <- log(y[i]);
    log_1my[i] <- log1m(y[i]);
  }
}
parameters {
  vector[K] gamma;
  real<lower=0> phi;
}
model {
  vector[N] Xgamma;
  real mu;
  real alpha_m1;
  real beta_m1;
  Xgamma <- X * gamma;
  for (i in 1:N) {
    mu <- inv_logit(Xgamma[i]);
    alpha_m1 <- mu * phi - 1.0;
    beta_m1 <- (1.0 - mu) * phi - 1.0;
    lp__ <- lp__ - lbeta(alpha,beta) + alpha_m1 * log_y[i] + 
                                        beta_m1 * log_1my[i];
  }
  // optional priors on gamma and phi here
}