Régression généralisée pondérée dans BUGS, JAGS

10

Dans Ron peut "pondérer préalablement" une glmrégression via le paramètre des poids . Par exemple:

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

Comment cela peut-il être accompli dans un modèle JAGSou BUGS?

J'ai trouvé un article sur ce sujet, mais aucun ne fournit d'exemple. Je m'intéresse principalement aux exemples de Poisson et de régression logistique.

user28937
la source
+1 très bonne question! Je posais cette question à certains spécialistes bayésiens et ils disent seulement que dans certains cas (poids selon la covariable catégorielle), vous pouvez calculer la distribution postérieure des paramètres pour chaque catégorie, puis les combiner dans une moyenne pondérée. Ils ne m'ont cependant pas donné de solution générale, donc je serais vraiment intéressé si elle existe ou non!
Curieux

Réponses:

7

Il pourrait être tard ... mais,

Veuillez noter 2 choses:

  • L'ajout de points de données n'est pas conseillé car cela modifierait les degrés de liberté. Les estimations moyennes de l'effet fixe pourraient être bien estimées, mais toute inférence devrait être évitée avec de tels modèles. Il est difficile de «laisser parler les données» si vous les modifiez.
  • Bien sûr, cela ne fonctionne qu'avec des poids à valeur entière (vous ne pouvez pas dupliquer 0,5 point de données), ce qui n'est pas ce qui se fait dans la régression la plus pondérée (lm). En général, une pesée est créée sur la base de la variabilité locale estimée à partir de répétitions (par exemple 1 / s ou 1 / s ^ 2 à un «x» donné) ou sur la base de la hauteur de réponse (par exemple 1 / Y ou 1 / Y ^ 2, à un «x» donné).

Dans Jags, Bugs, Stan, proc MCMC, ou en bayésien en général, la probabilité n'est pas différente que dans lm ou glm fréquentiste (ou n'importe quel modèle), c'est tout de même !! Créez simplement une nouvelle colonne "poids" pour votre réponse et écrivez la probabilité

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

Ou un poisson lesté:

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

Ce code Bugs / Jags ferait tout simplement l'affaire. Vous obtiendrez tout correct. N'oubliez pas de continuer à multiplier le postérieur de tau par le poids, par exemple lors de la prédiction et des intervalles de confiance / prédiction.

Pierre Lebrun
la source
Si nous le faisons comme indiqué, nous modifions la moyenne et la variance. Pourquoi n'est-ce pas y [i] * poids [i] ~ dpois (lambda [i] * poids [i])? Cela ne ferait que régler la variance. Le problème ici est que y [i] * weight [i] peut être de type réel.
user28937
en effet, la régression pondérée change de moyen (car la pesée conduit la régression à se rapprocher des points qui ont beaucoup de poids!) et la variance est désormais fonction des poids (donc ce n'est pas un modèle homoskédastique). La variance (ou précision) tau n'a plus de sens, mais tau / poids [i] peut être interprété exactement comme la précision du modèle (pour un "x" donné). Je ne conseillerais pas la multiplication des données (y) par les poids ... attendez-vous si c'est précisément quelque chose que vous voulez faire, mais je ne comprends pas votre modèle dans ce cas ...
Pierre Lebrun
Je suis d'accord avec vous, cela ne change pas la moyenne dans l'exemple normal: y [i] ~ dnorm (mu [i], tau / poids [i]), mais c'est le cas dans le second, puisque lambda [i] * poids [ i] devient la "nouvelle" lambda pour dpois et cela ne correspondra plus à y [i]. Je dois me corriger, cela devrait être: ty [i] * exp (poids [i]) ~ dpois (lambda [i] * poids [i]). L'idée avec la multiplication dans le cas du poisson est que nous voulons ajuster la variance, mais aussi ajuster la moyenne, donc n'avons-nous pas besoin de corriger la moyenne?
user28937
Si vous devez ajuster la variance indépendamment, peut-être qu'un modèle binomial négatif peut être utile au lieu d'un Poisson? Il ajoute un paramètre d'inflation / déflation de la variance au Poisson. Sauf que c'est très similaire.
Pierre Lebrun
Pierre bonne idée. J'ai aussi pensé à la représentation continue de la distribution de poisson définie dans la diapositive 6/12 sur linkd
user28937
4

Tout d'abord, il convient de souligner qu'il glmn'effectue pas de régression bayésienne. Le paramètre «poids» est essentiellement une abréviation de «proportion d'observations», qui peut être remplacée par un suréchantillonnage de votre ensemble de données de manière appropriée. Par exemple:

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)

Donc, pour ajouter du poids aux points dans JAGS ou BUGS, vous pouvez augmenter votre jeu de données de la même manière que ci-dessus.

David Marx
la source
2
cela ne fonctionnera pas, les poids sont généralement des nombres réels, pas des nombres entiers
Curieux
Cela ne vous empêche pas de les approximer avec des nombres entiers. Ma solution n'est pas parfaite, mais elle fonctionne approximativement. Par exemple, compte tenu des poids (1/3, 2/3, 1), vous pouvez suréchantillonner la deuxième classe d'un facteur deux et la troisième classe d'un facteur trois.
David Marx
0

J'ai essayé d'ajouter un commentaire ci-dessus, mais mon représentant est trop faible.

Devrait

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

ne pas être

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

dans JAGS? J'exécute des tests comparant les résultats de cette méthode dans JAGS aux résultats d'une régression pondérée via lm () et ne peut trouver de conformité qu'en utilisant cette dernière. Voici un exemple simple:

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)

et comparer à

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
la source
Indépendamment de la réputation, les commentaires ne doivent pas être donnés comme réponses.
Michael R. Chernick