Как я могу смоделировать пропорции с BUGS / JAGS / STAN?

10

Я пытаюсь построить модель, в которой ответ является пропорцией (фактически это доля голосов, которую партия получает в избирательных округах). Его распределение не является нормальным, поэтому я решил смоделировать его с помощью бета-версии. У меня также есть несколько предикторов.

Однако я не знаю, как написать это в BUGS / JAGS / STAN (JAGS был бы моим лучшим выбором, но это не имеет значения). Моя проблема в том, что я делаю сумму параметров по предикторам, но что я могу с этим сделать?

Код будет что - то вроде этого (в синтаксисе зазубрины), но знает , что Дон»как„звено“ y_hatи yпараметры.

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

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

( y_hatэто просто перекрестное произведение параметров и предикторов, отсюда и детерминированные отношения. aИ bэто коэффициенты, которые я пытаюсь оценить, xбудучи предиктором).

Спасибо за ваши предложения!

Жоэль
источник
Что такое a, b, y_hat? Вы должны четко определить свою модель. Кстати, синтаксис BUGS близок к математическому синтаксису. Таким образом, если вы знаете, как написать свою модель на математическом языке, то почти вся работа выполнена.
Стефан Лоран
Стефан, спасибо. Я отредактировал вопрос, чтобы определить a, b, y_hat. Я также не знаю математического ответа, иначе ответ был бы намного проще ;-)
Joël
Я подозреваю, что я мог бы опираться на тот факт, что E (y) = альфа / (альфа + бета), но я не могу понять, как именно.
Жоэль

Ответы:

19

Подход бета-регрессии заключается в повторной параметризации в терминах и . Где будет эквивалентно y_hat, которое вы предсказываете. В этой параметризации у вас будет и . Затем вы можете смоделировать как логит линейной комбинации. может иметь собственный априор (должен быть больше 0) или может быть смоделирован также на ковариатах (выберите функцию связи, чтобы оставить ее больше 0, например, экспоненциальную).μφμαзнак равноμ×φβзнак равно(1-μ)×φμφ

Возможно что-то вроде:

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)
Грег Сноу
источник
Спасибо, это очень полезно! Я пытаюсь подобрать модель по твоему совету.
Жоэль
Однако, когда я запускаю модель, я получаю ошибки, такие как: «Ошибка в узле y [6283] Недопустимые родительские значения». Есть идеи, что здесь происходит?
Жоэль
@ Joël, какое значение у [6283]? Вы убедились, что значения альфа и бета-версий ограничены допустимыми значениями? Я ожидаю, что что-то может пойти до 0 или ниже, и это дает ошибку.
Грег Сноу
Нет, я проверил это, все мои значения y строго превосходят 0 (и уступают 1). Может быть, мои приоры сталкиваются с эмпирическими значениями y в какой-то момент? Но я не знаю, как это проверить, и мои приоры кажутся разумными - по крайней мере, мне!
Жоэль
1
@colin, я плохо знаю JAGS, так что лучше спросить об этом на форуме специально для JAGS. Или попробуйте другой инструмент, я считаю, что в эти дни мне нравится Stan for Bayes.
Грег Сноу
18

Грег Сноу дал отличный ответ. Для полноты здесь приведен эквивалент в синтаксисе Stan. Хотя у Stan есть бета-дистрибутив, который вы можете использовать, быстрее вычислить логарифм бета-плотности самостоятельно, потому что константы log(y)и log(1-y)могут быть рассчитаны один раз в начале (а не каждый раз, когда y ~ beta(alpha,beta)это будет вызываться). Увеличивая зарезервированную lp__переменную (см. Ниже), вы можете суммировать логарифм бета-плотности по наблюдениям в вашей выборке. Я использую метку «гамма» для вектора параметров в линейном предикторе.

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
}
Бен Гудрич
источник
Спасибо Бен! Очень полезно иметь синтаксис Stan.
Жоэль
У Стэна v2 есть оператор выборки "beta_proportion", который, я считаю, устраняет необходимость непосредственно манипулировать "lp__"
THK