BUGS、JAGS 中的加权广义回归

机器算法验证 贝叶斯 广义线性模型 锯齿 错误 加权回归
2022-03-04 04:18:46

我们可以通过weights参数R“先验权重”glm回归。例如:

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

这如何在 a JAGSorBUGS模型中完成?

我找到了一些讨论这个问题的论文,但没有一个提供例子。我主要对泊松和逻辑回归示例感兴趣。

3个回答

可能会迟到……但是,

请注意两点:

  • 不建议添加数据点,因为它会改变自由度。可以很好地估计固定效应的平均估计值,但应避免使用此类模型进行所有推断。如果你改变它,就很难“让数据说话”。
  • 当然,它仅适用于整数值权重(您不能复制 0.5 个数据点),这不是大多数加权 (lm) 回归中所做的。一般来说,权重是基于从重复估计的局部变异性(例如,给定“x”处的 1/s 或 1/s^2)或基于响应高度(例如,1/Y 或 1/Y^2,在给定的'x')。

在 Jags、Bugs、Stan、proc MCMC 或一般的贝叶斯中,可能性与常客 lm 或 glm(或任何模型)没有什么不同,它是一样的!只需为您的响应创建一个新列“权重”,并将可能性写为

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

或加权泊松:

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

这个 Bugs/Jags 代码很简单。你会得到一切正确的。不要忘记继续将 tau 的后验乘以权重,例如在进行预测和置信/预测区间时。

首先,值得指出的是glm不执行贝叶斯回归。“权重”参数基本上是“观察比例”的简写,可以用适当地对数据集进行上采样来代替。例如:

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)

因此,要为 JAGS 或 BUGS 中的点添加权重,您可以以与上述类似的方式扩充您的数据集。

尝试在上面添加评论,但我的代表太低了。

应该

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

不是

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

在JAGS?我正在运行一些测试,将 JAGS 中这种方法的结果与通过 lm() 加权回归的结果进行比较,并且只能使用后者找到一致性。这是一个简单的例子:

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)

并比较

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