Regresión generalizada ponderada en BUGS, JAGS

En Rpodemos "ponderar previamente" una glmregresión a través del parámetro de pesos . Por ejemplo:

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

¿Cómo se puede lograr esto en una JAGSo BUGSmodelo?

Encontré un artículo discutiendo esto, pero ninguno de ellos proporciona un ejemplo. Me interesan principalmente los ejemplos de regresión logística y Poisson.

Puede que sea tarde ... pero,

Tenga en cuenta 2 cosas:

• No se recomienda agregar puntos de datos, ya que cambiaría los grados de libertad. Las estimaciones medias de efecto fijo podrían estimarse bien, pero se debe evitar toda inferencia con tales modelos. Es difícil "dejar que los datos hablen" si los cambia.
• Por supuesto, solo funciona con pesos con valores enteros (no puede duplicar 0,5 puntos de datos), que no es lo que se hace en la regresión más ponderada (lm). En general, se crea un pesaje basado en la variabilidad local estimada a partir de réplicas (por ejemplo, 1 / so 1 / s ^ 2 en una 'x' dada) o en función de la altura de respuesta (por ejemplo, 1 / Y o 1 / Y ^ 2, en una 'x' dada).

En Jags, Bugs, Stan, proc MCMC, o en Bayesian en general, la probabilidad no es diferente que en la frecuente lm o glm (o cualquier modelo), ¡¡¡es lo mismo !! Simplemente cree una nueva columna "ponderación" para su respuesta y escriba la probabilidad como

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

O un poisson ponderado:

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

Este código Bugs / Jags sería simplemente el truco. Obtendrás todo correcto. No olvide continuar multiplicando la parte posterior de tau por el peso, por ejemplo, al hacer predicciones e intervalos de confianza / predicción.

Primero, vale la pena señalar que glmno realiza una regresión bayesiana. El parámetro 'pesos' es básicamente una abreviatura de "proporción de observaciones", que se puede reemplazar con un muestreo ascendente de su conjunto de datos de manera adecuada. Por ejemplo:

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)

Entonces, para agregar peso a los puntos en JAGS o BUGS, puede aumentar su conjunto de datos de manera similar a la anterior.

Intenté agregar al comentario anterior, pero mi representante es demasiado bajo.

Debería

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

No ser

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

en JAGS? Estoy ejecutando algunas pruebas que comparan los resultados de este método en JAGS con los resultados de una regresión ponderada a través de lm () y solo puedo encontrar la conformidad con este último. Aquí hay un ejemplo 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)

y comparar con

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