Gewichtete generalisierte Regression in BUGS, JAGS

10

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.

user28937
quelle
+1 sehr gute Frage! Ich habe dies einige Bayes'sche Spezialisten gefragt und sie sagen nur, dass Sie in einigen Fällen (Gewichte nach kategorialer Kovariate) die posteriore Verteilung der Parameter für jede Kategorie berechnen und sie dann in einem gewichteten Mittelwert kombinieren können. Sie gaben mir jedoch keine allgemeine Lösung, daher wäre ich wirklich interessiert, ob es existiert oder nicht!
Neugierig

Antworten:

7

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.

Pierre Lebrun
quelle
Wenn wir es wie angegeben tun, ändern wir den Mittelwert und die Varianz. Warum ist es nicht y [i] * weight [i] ~ dpois (Lambda [i] * weight [i])? Das würde nur die Varianz anpassen. Das Problem hierbei ist, dass y [i] * weight [i] vom Typ real sein kann.
user28937
In der Tat ändert sich die gewichtete Regression im Mittel (weil das Wiegen dazu führt, dass die Regression näher an die Punkte mit vielen Gewichten heranrückt!) und die Varianz nun eine Funktion der Gewichte ist (daher handelt es sich nicht um ein homoskedastisches Modell). Die Varianz (oder Präzision) Tau hat keine Bedeutung mehr, aber Tau / Gewicht [i] kann genau als Präzision des Modells interpretiert werden (für ein gegebenes "x"). Ich würde die Multiplikation der Daten (y) mit den Gewichten nicht empfehlen ... erwarten Sie, ob dies genau etwas ist, was Sie tun möchten, aber ich verstehe Ihr Modell in diesem Fall nicht ...
Pierre Lebrun
Ich stimme Ihnen zu, dass es den Mittelwert im normalen Beispiel nicht ändert: y [i] ~ dnorm (mu [i], tau / weight [i]), aber im zweiten, da Lambda [i] * weight [ i] wird das "neue" Lambda für dpois und dies wird nicht mehr mit y [i] übereinstimmen. Ich muss mich korrigieren, es sollte sein: ty [i] * exp (Gewicht [i]) ~ dpois (Lambda [i] * Gewicht [i]). Die Idee mit der Multiplikation im Poisson-Fall ist, dass wir die Varianz anpassen wollen, aber auch den Mittelwert, also müssen wir den Mittelwert nicht korrigieren?
user28937
Wenn Sie die Varianz unabhängig anpassen müssen, kann möglicherweise ein negatives Binomialmodell anstelle eines Poisson nützlich sein? Es fügt dem Poisson einen Varianzinflations- / Deflationsparameter hinzu. Nur dass es sehr ähnlich ist.
Pierre Lebrun
Pierre gute Idee. Ich dachte auch an die kontinuierliche Darstellung der Poisson-Verteilung, die in Folie 6/12 unter linkd
user28937 am
4

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.

David Marx
quelle
2
Dies wird nicht funktionieren, Gewichte sind normalerweise reelle Zahlen, keine ganzen Zahlen
Neugierig
Dies schließt nicht aus, dass Sie sie mit ganzen Zahlen approximieren. Meine Lösung ist nicht perfekt, funktioniert aber ungefähr. Bei Gewichten (1/3, 2/3, 1) können Sie beispielsweise die zweite Klasse um den Faktor zwei und die dritte Klasse um den Faktor drei erhöhen.
David Marx
0

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
metrikrn
quelle
Unabhängig von der Reputation sollten Kommentare nicht als Antworten gegeben werden.
Michael R. Chernick