Ważona regresja uogólniona w BŁĘDACH, JAGACH

10

W R„uprzedniej wadze” możemy glmregresję za pomocą parametru wag . Na przykład:

glm.D93 <- glm(counts ~ outcome + treatment, family = poisson(), weights=w)

Jak można tego dokonać w modelu JAGSlub BUGS?

Znalazłem jakiś artykuł na ten temat, ale żaden z nich nie stanowi przykładu. Interesują mnie głównie przykłady Poissona i regresji logistycznej.

użytkownik 28937
źródło
+1 bardzo dobre pytanie! Pytałem o to kilku bayesowskich specjalistów, którzy mówią tylko, że w niektórych przypadkach (wagi według współzmienności kategorialnej) można obliczyć rozkład tylny parametrów dla każdej kategorii, a następnie połączyć je w średnią ważoną. Nie dali mi jednak ogólnego rozwiązania, więc byłbym bardzo zainteresowany, czy istnieje, czy nie!
Ciekawy

Odpowiedzi:

7

Może być późno ... ale ...

Uwaga 2 rzeczy:

  • Dodawanie punktów danych nie jest zalecane, ponieważ zmieniłoby to stopnie swobody. Średnie oszacowania ustalonego efektu można dobrze oszacować, ale takich modeli należy unikać. Trudno jest „pozwolić, aby dane mówiły”, jeśli je zmienisz.
  • Oczywiście działa tylko z wagami o wartościach całkowitych (nie można powielić 0,5 punktu danych), co nie jest wykonywane w większości regresji ważonej (lm). Zasadniczo ważenie jest tworzone na podstawie lokalnej zmienności oszacowanej na podstawie powtórzeń (np. 1 / s lub 1 / s ^ 2 przy danym „x”) lub na podstawie wysokości odpowiedzi (np. 1 / Y lub 1 / Y ^ 2, przy dany „x”).

W Jags, Bugs, Stan, proc MCMC lub w Bayesian w ogóle prawdopodobieństwo nie jest inne niż w lmistom lm ​​lub glm (lub dowolnym innym modelu), jest po prostu takie samo !! Po prostu utwórz nową kolumnę „waga” dla swojej odpowiedzi i wpisz prawdopodobieństwo jako

y [i] ~ dnorm (mu [i], tau / weight [i])

Lub ważona poissona:

y [i] ~ dpois (lambda [i] * waga [i])

Ten kod Bugs / Jags po prostu przydałby się. Dostaniesz wszystko poprawnie. Nie zapomnij o dalszym pomnożeniu tylnej części tau przez wagę, na przykład przy dokonywaniu prognoz i przedziałach ufności / prognozy.

Pierre Lebrun
źródło
Jeśli zrobimy to, jak powiedziano, zmieniamy średnią i wariancję. Dlaczego to nie y [i] * waga [i] ~ dpois (lambda [i] * waga [i])? To skorygowałoby tylko wariancję. Problem polega na tym, że y [i] * waga [i] może być typu rzeczywistego.
user28937
w rzeczywistości regresja ważona zmienia się jako średnia (ponieważ ważenie prowadzi regresję do przybliżenia punktów, które mają wiele wag!), a wariancja jest teraz funkcją wag (stąd nie jest to model homoskedastyczny). Wariancja (lub precyzja) tau nie ma już znaczenia, ale tau / waga [i] może być interpretowana dokładnie jako dokładność modelu (dla danego „x”). Nie radziłbym mnożenia danych (y) przez wagi ... spodziewaj się, że jest to dokładnie to, co chcesz zrobić, ale nie rozumiem twojego modelu w tym przypadku ...
Pierre Lebrun
Zgadzam się z tobą, że nie zmienia to średniej w normalnym przykładzie: y [i] ~ dnorm (mu [i], tau / weight [i]), ale robi to w drugim, ponieważ lambda [i] * waga [ i] staje się „nową” lambda dla dpois i to już nie będzie pasować do y [i]. Muszę się poprawić to powinno być: ty [i] * exp (waga [i]) ~ dpois (lambda [i] * waga [i]). Pomysł z pomnożeniem w przypadku Poissona polega na tym, że chcemy dostosować wariancję, ale także dostosować średnią, więc czy nie musimy korygować średniej?
user28937
Jeśli musisz samodzielnie dopasować wariancję, może przydatny będzie model dwumianowy ujemny zamiast Poissona? Dodaje parametr inflacji / deflacji wariancji do Poissona. Tyle że jest bardzo podobny.
Pierre Lebrun
Pierre dobry pomysł. Myślałem również o ciągłej reprezentacji rozkładu Poissona
28937
4

Po pierwsze, warto zauważyć, że glmnie wykonuje regresji bayesowskiej. Parametr „wagi” to w skrócie „proporcja obserwacji”, którą można zastąpić odpowiednim próbkowaniem zestawu danych. Na przykład:

x=1:10
y=jitter(10*x)
w=sample(x,10)

augmented.x=NULL
augmented.y=NULL    
for(i in 1:length(x)){
    augmented.x=c(augmented.x, rep(x[i],w[i]))
    augmented.y=c(augmented.y, rep(y[i],w[i]))
}

# These are both basically the same thing
m.1=lm(y~x, weights=w)
m.2=lm(augmented.y~augmented.x)

Aby dodać wagę do punktów w JAGS lub BŁĘDACH, możesz rozszerzyć swój zestaw danych w podobny sposób jak powyżej.

David Marks
źródło
2
to nie zadziała, wagi są zwykle liczbami rzeczywistymi, a nie liczbami całkowitymi
Ciekawy
Nie wyklucza to przybliżenia ich liczbami całkowitymi. Moje rozwiązanie nie jest idealne, ale działa w przybliżeniu. Na przykład, biorąc pod uwagę wagi (1/3, 2/3, 1), możesz upsamplować drugą klasę dwa razy, a trzecią klasę trzy razy.
David Marx,
0

Próbowałem dodać do komentarza powyżej, ale moje powtórzenie jest zbyt niskie.

Powinien

y[i] ~ dnorm(mu[i], tau / weight[i])

nie być

y[i] ~ dnorm(mu[i], tau * weight[i])

w JAGS? Przeprowadzam kilka testów porównujących wyniki tej metody w JAGS z wynikami ważonej regresji za pomocą lm () i mogę znaleźć zgodność tylko przy użyciu tej ostatniej. Oto prosty przykład:

aggregated <- 
  data.frame(x=1:5) %>%
  mutate( y = round(2 * x + 2 + rnorm(length(x)) ),
          freq = as.numeric(table(sample(1:5, 100, 
                 replace=TRUE, prob=c(.3, .4, .5, .4, .3)))))
x <- aggregated$x
y <- aggregated$y
weight <- aggregated$freq
N <- length(y)

# via lm()
lm(y ~ x, data = aggregated, weight = freq)

i porównaj z

lin_wt_mod <- function() {

  for (i in 1:N) {
    y[i] ~ dnorm(mu[i], tau*weight[i])
    mu[i] <- beta[1] + beta[2] * x[i]
  }

  for(j in 1:2){
    beta[j] ~ dnorm(0,0.0001)
  }

  tau   ~ dgamma(0.001, 0.001)
  sigma     <- 1/sqrt(tau)
}

dat <- list("N","x","y","weight")
params <- c("beta","tau","sigma")

library(R2jags)
fit_wt_lm1 <- jags.parallel(data = dat, parameters.to.save = params,
              model.file = lin_wt_mod, n.iter = 3000, n.burnin = 1000)
fit_wt_lm1$BUGSoutput$summary
metrikrn
źródło
Niezależnie od reputacji komentarze nie powinny być podawane jako odpowiedzi.
Michael R. Chernick