Jak mogę modelować proporcje za pomocą BŁĘDÓW / JAGS / STAN?

Próbuję zbudować model, w którym odpowiedź jest proporcjonalna (w rzeczywistości jest to liczba głosów, jaką partia zdobywa w okręgach wyborczych). Jego rozkład nie jest normalny, więc postanowiłem modelować go z rozkładem beta. Mam też kilka predyktorów.

Nie wiem jednak, jak napisać to w BŁĘDACH / JAGACH / STANU (JAGS byłby moim najlepszym wyborem, ale to naprawdę nie ma znaczenia). Mój problem polega na tym, że wykonuję sumę parametrów za pomocą predyktorów, ale co mogę z tym zrobić?

Kod byłby mniej więcej taki (w składni JAGS), ale nie wiem, jak „połączyć” parametry y_hati y.

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

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

( y_hatjest tylko produktem krzyżowym parametrów i predyktorów, stąd relacja deterministyczna. ai bsą współczynnikami, które staram się oszacować, xbędąc predyktorem).

Dzięki za sugestie!

— Joël
źródło

Co to jest a, b, y_hat? Powinieneś jasno zdefiniować swój model. Nawiasem mówiąc, składnia BŁĘDÓW jest zbliżona do składni matematycznej. Zatem jeśli wiesz, jak napisać swój model w języku matematycznym, prawie cała praca zostanie wykonana.

— Stéphane Laurent,

Stéphane, dzięki. Zredagowałem pytanie, aby zdefiniować a, b, y_hat. Matematycznie też nie znam odpowiedzi, w przeciwnym razie odpowiedź byłaby znacznie łatwiejsza ;-)

— Joël

Podejrzewam, że mógłbym oprzeć się na fakcie, że E (y) = alpha / (alpha + beta), ale tak naprawdę nie jestem w stanie dowiedzieć się, jak dokładnie.

— Joël

Odpowiedzi:

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

Być może coś takiego:

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
źródło

Dziękuję, to jest bardzo pomocne! Staram się dopasować model do twojej rady.

— Joël

Jednak po uruchomieniu modelu pojawiają się błędy, takie jak: „Błąd w węźle y [6283] Nieprawidłowe wartości nadrzędne”. Masz pomysł, co się tutaj dzieje?

— Joël

@ Joël, jaka jest wartość y [6283]? czy upewniłeś się, że wartości alf i bet są ograniczone do wartości prawnych? Spodziewam się, że coś mogło spaść do 0 lub poniżej i to powoduje błąd.

— Greg Snow,

Nie, sprawdziłem, że wszystkie moje wartości y są ściśle wyższe od 0 (i niższe od 1). Może w pewnym momencie moje priory ścierają się z empirycznymi wartościami y? Ale nie wiem, jak to sprawdzić, a moje przeory wydają się rozsądne - przynajmniej dla mnie!

— Joël

@colin, nie znam tak dobrze JAGS, więc lepiej zadać to pytanie na forum specjalnie dla JAGS. Lub wypróbuj to w innym narzędziu, obecnie lubię Stana za Bayesa.

— Greg Snow,

Greg Snow dał świetną odpowiedź. Dla kompletności, oto odpowiednik w składni Stan. Chociaż Stan ma rozkład beta, którego można użyć, szybciej jest samodzielnie opracować logarytm gęstości beta, ponieważ stałe log(y)i log(1-y)można je obliczyć na początku (a nie za każdym razem, gdy y ~ beta(alpha,beta)zostanie to wywołane). Zwiększając zarezerwowaną lp__zmienną (patrz poniżej), można zsumować logarytm gęstości beta z obserwacji w próbie. Używam etykiety „gamma” dla wektora parametru w predyktorze liniowym.

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
źródło

Dzięki Ben! Bardzo przydatne, aby mieć również składnię Stan.

— Joël

Stan v2 ma oświadczenie próbkowania „beta_proporcja”, które moim zdaniem eliminuje potrzebę bezpośredniego manipulowania „lp__”

— THK