10
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ What are a, b, y_hat ? You should clearly define your model. By the way the BUGS syntax is close to the mathematical syntax. Thus if you know how to write your model in mathematical language then almost all the job is done. $\endgroup$ – Stéphane Laurent Oct 30 '12 at 19:27
  • $\begingroup$ Stéphane, thanks. I edited the question to define a, b, y_hat. I don't know the answer mathematically either, otherwise the answer would indeed be much easier ;-) $\endgroup$ – Joël Oct 30 '12 at 19:33
  • $\begingroup$ I suspect that I could build on the fact that E(y) = alpha / (alpha + beta), but I can't really figure out how exactly. $\endgroup$ – Joël Oct 30 '12 at 19:35
19
$\begingroup$

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)
$\endgroup$
  • $\begingroup$ Thank you, this is very helpful! I'm trying to fit a model with your advice. $\endgroup$ – Joël Oct 30 '12 at 20:45
  • $\begingroup$ However, when I run the model, I get errors such as : "Error in node y[6283] Invalid parent values". Any idea what is going on here? $\endgroup$ – Joël Oct 30 '12 at 20:46
  • $\begingroup$ @Joël, what is the value of y[6283]? have you made sure that the values of the alphas and betas are restricted to legal values? I expect that something may have gone to 0 or below and that gives the error. $\endgroup$ – Greg Snow Oct 31 '12 at 15:04
  • $\begingroup$ No, I checked that, all my y values are strictly superior to 0 (and inferior to 1). Maybe my priors clash with the empirical y values at some point? But I don't know how to check this, and my priors seem sensible - at least to me! $\endgroup$ – Joël Oct 31 '12 at 15:16
  • 1
    $\begingroup$ @colin, I don't know JAGS that well, so this may be better asked on a forum specifically for JAGS. Or try it in a different tool, I find that I like Stan for Bayes these days. $\endgroup$ – Greg Snow Aug 17 '16 at 15:52
18
$\begingroup$

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
}
$\endgroup$
  • $\begingroup$ Thanks Ben! Very useful to have the Stan syntax as well. $\endgroup$ – Joël Nov 6 '12 at 23:59
  • $\begingroup$ Stan v2 has a "beta_proportion" sampling statement that I believe obviates the need to directly manipulate "lp__" $\endgroup$ – THK Jun 2 '19 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.