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

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.

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.

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.

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