Como posso modelar uma proporção com BUGS / JAGS / STAN?

Estou tentando construir um modelo em que a resposta seja proporcional (na verdade, é a parcela de votos que um partido obtém nos círculos eleitorais). Sua distribuição não é normal, então decidi modelá-la com uma distribuição beta. Eu também tenho vários preditores.

No entanto, não sei como escrevê-lo em BUGS / JAGS / STAN (JAGS seria minha melhor escolha, mas isso realmente não importa). Meu problema é que eu faço uma soma de parâmetros pelos preditores, mas o que posso fazer com isso?

O código seria algo assim (na sintaxe JAGS), mas não sei como "vincular" os parâmetros y_hate y.

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

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

( y_haté apenas o produto cruzado de parâmetros e preditores, daí a relação determinística. ae bsão os coeficientes que tento estimar, xsendo preditores).

Obrigado por suas sugestões!

$\mu$$\phi$$\mu$$\alpha=\mu\times\phi$$\beta=(1-\mu) \times \phi$$\mu$$\phi$

Possivelmente algo como:

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)

Greg Snow deu uma ótima resposta. Para completar, aqui está o equivalente na sintaxe Stan. Embora Stan tenha uma distribuição beta que você possa usar, é mais rápido calcular o logaritmo da densidade beta, porque as constantes log(y)e log(1-y)podem ser calculadas uma vez desde o início (e não sempre que y ~ beta(alpha,beta)seria chamado). Ao incrementar a lp__variável reservada (veja abaixo), você pode somar o logaritmo da densidade beta sobre as observações em sua amostra. Eu uso o rótulo "gama" para o vetor de parâmetro no preditor linear.

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
}