Gewichtete generalisierte Regression in BUGS, JAGS

In können Rwir eine glmRegression über den Gewichtungsparameter "vorgewichten" . Zum Beispiel:

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

Wie kann dies in einem JAGSoder BUGSModell erreicht werden?

Ich habe ein Papier gefunden, das dies diskutiert, aber keines davon liefert ein Beispiel. Ich interessiere mich hauptsächlich für Poisson und logistische Regressionsbeispiele.

Es könnte spät sein ... aber,

Bitte beachten Sie 2 Dinge:

• Das Hinzufügen von Datenpunkten wird nicht empfohlen, da dies die Freiheitsgrade ändern würde. Mittlere Schätzungen des festen Effekts könnten gut geschätzt werden, aber alle Schlussfolgerungen sollten mit solchen Modellen vermieden werden. Es ist schwierig, "die Daten sprechen zu lassen", wenn Sie sie ändern.
• Natürlich funktioniert es nur mit ganzzahligen Gewichten (Sie können 0,5 Datenpunkte nicht duplizieren), was bei der am meisten gewichteten (lm) Regression nicht der Fall ist. Im Allgemeinen wird eine Wägung basierend auf der lokalen Variabilität erstellt, die aus Wiederholungen (z. B. 1 / s oder 1 / s ^ 2 bei einem gegebenen 'x') geschätzt wird, oder basierend auf der Antworthöhe (z. B. 1 / Y oder 1 / Y ^ 2 bei) ein gegebenes 'x').

In Jags, Bugs, Stan, Proc MCMC oder in Bayesian im Allgemeinen ist die Wahrscheinlichkeit nicht anders als in Frequentist lm oder glm (oder einem anderen Modell), es ist genau das gleiche !! Erstellen Sie einfach eine neue Spalte "Gewicht" für Ihre Antwort und schreiben Sie die Wahrscheinlichkeit als

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

Oder eine gewichtete Poisson:

y [i] ~ dpois (Lambda [i] * Gewicht [i])

Dieser Bugs / Jags-Code würde einfach zum Trick. Sie werden alles richtig machen. Vergessen Sie nicht, den hinteren Teil von Tau weiter mit dem Gewicht zu multiplizieren, beispielsweise wenn Sie Vorhersagen und Konfidenz- / Vorhersageintervalle treffen.

Zunächst ist darauf hinzuweisen, dass glmkeine Bayes'sche Regression durchgeführt wird. Der Parameter 'Gewichte' ist im Grunde eine Abkürzung für "Anteil der Beobachtungen", die durch eine entsprechende Aktualisierung Ihres Datensatzes ersetzt werden kann. Zum Beispiel:

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)

Um Punkte in JAGS oder BUGS zu gewichten, können Sie Ihren Datensatz auf ähnliche Weise wie oben erweitern.

Ich habe versucht, den obigen Kommentar hinzuzufügen, aber mein Repräsentant ist zu niedrig.

Sollte

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

nicht sein

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

in JAGS? Ich führe einige Tests durch, in denen die Ergebnisse dieser Methode in JAGS mit den Ergebnissen einer gewichteten Regression über lm () verglichen werden, und kann nur mit letzterer Übereinstimmung finden. Hier ist ein einfaches Beispiel:

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)

und vergleiche mit

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