¿Cómo puedo modelar una proporción con BUGS / JAGS / STAN?

10

Estoy tratando de construir un modelo donde la respuesta sea proporcional (en realidad es la proporción de votos que un partido obtiene en los distritos electorales). Su distribución no es normal, así que decidí modelarla con una distribución beta. También tengo varios predictores.

Sin embargo, no sé cómo escribirlo en BUGS / JAGS / STAN (JAGS sería mi mejor opción, pero en realidad no importa). Mi problema es que hago una suma de parámetros por predictores, pero ¿qué puedo hacer con ellos?

El código sería algo como esto (en la sintaxis JAGS), pero no sé cómo "vincular" los parámetros y_haty y.

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

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

( y_hates solo el producto cruzado de parámetros y predictores, de ahí la relación determinista. ay bson los coeficientes que trato de estimar, xsiendo un predictor).

Gracias por tus sugerencias!

Joël
fuente
¿Qué son a, b, y_hat? Debe definir claramente su modelo. Por cierto, la sintaxis de BUGS está cerca de la sintaxis matemática. Por lo tanto, si sabe cómo escribir su modelo en lenguaje matemático, casi todo el trabajo está hecho.
Stéphane Laurent
Stéphane, gracias. Edité la pregunta para definir a, b, y_hat. Tampoco sé la respuesta matemáticamente, de lo contrario la respuesta sería mucho más fácil ;-)
Joël
Sospecho que podría construir sobre el hecho de que E (y) = alfa / (alpha + beta), pero realmente no puedo entender cómo exactamente.
Joël

Respuestas:

19

El enfoque de regresión beta es volver a parametrizar en términos de y . Donde será el equivalente a y_hat que predices. En esta parametrización tendrá y . Luego puede modelar como el logit de la combinación lineal. puede tener su propio previo (debe ser mayor que 0), o también puede modelarse en covariables (elija una función de enlace para mantenerlo mayor que 0, como exponencial).ϕ μ α = μ × ϕ β = ( 1 - μ ) × ϕ μ ϕμϕμα=μ×ϕβ=(1μ)×ϕμϕ

Posiblemente 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
fuente
¡Gracias, esto es muy útil! Estoy tratando de adaptar un modelo con tu consejo.
Joël
Sin embargo, cuando ejecuto el modelo, obtengo errores como: "Error en el nodo y [6283] Valores primarios no válidos". ¿Alguna idea de lo que está pasando aquí?
Joël
@ Joël, ¿cuál es el valor de y [6283]? ¿Se ha asegurado de que los valores de los alfa y las beta estén restringidos a valores legales? Espero que algo haya ido a 0 o menos y eso da el error.
Greg Snow
No, lo comprobé, todos mis valores y son estrictamente superiores a 0 (e inferiores a 1). ¿Quizás mis antecedentes chocan con los valores empíricos y en algún momento? Pero no sé cómo verificar esto, y mis antecedentes parecen ser sensatos, ¡al menos para mí!
Joël
1
@colin, no conozco JAGS tan bien, por lo que es mejor preguntar en un foro específicamente para JAGS. O pruébelo en una herramienta diferente, encuentro que me gusta Stan para Bayes en estos días.
Greg Snow
18

Greg Snow dio una gran respuesta. Para completar, aquí está el equivalente en la sintaxis de Stan. Aunque Stan tiene una distribución beta que podría usar, es más rápido calcular el logaritmo de la densidad beta usted mismo debido a las constantes log(y)y log(1-y)puede calcularse una vez desde el principio (en lugar de cada vez que y ~ beta(alpha,beta)se llamaría). Al incrementar la lp__variable reservada (ver más abajo), puede sumar el logaritmo de la densidad beta sobre las observaciones en su muestra. Utilizo la etiqueta "gamma" para el vector de parámetros en el predictor lineal.

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
}
Ben Goodrich
fuente
Gracias ben! Muy útil para tener la sintaxis de Stan también.
Joël
Stan v2 tiene una declaración de muestreo "beta_proportion" que creo que evita la necesidad de manipular directamente "lp__"
THK