拟合贝叶斯线性回归并预测不可观察的值

Posted 2023-02-19

【中文标题】拟合贝叶斯线性回归并预测不可观察的值

【英文标题】：Fit a bayesian linear regression and predict unobservable values

【发布时间】：2016-02-13 06:40:42

【问题描述】：

我想使用 Jags 加 R 来调整具有可观察量的线性模型，并对不可观察量进行推断。我在互联网上找到了很多关于如何调整模型的示例，但没有关于在 Jags 环境中拟合模型后如何推断其系数的示例。因此，我将不胜感激。

我的数据如下所示：

ngroups <- 2
group <- 1:ngroups
nobs <- 100
dta <- data.frame(group=rep(group,each=nobs),y=rnorm(nobs*ngroups),x=runif(nobs*ngroups))
head(dta)

【问题讨论】：

你真的想要一个信息丰富的先验吗？如果没有，只需使用lm 和predict 函数。

【参考方案1】：

JAGS 有强大的方法来推断丢失的数据，一旦掌握了它，就很容易了！我强烈建议您查看 Marc Kéry 的 excellent book，它提供了对 BUGS 语言编程的精彩介绍（JAGS 与 BUGS 非常接近，几乎所有内容都可以迁移）。

如您所说，最简单的方法是调整模型。下面我提供了一个完整的例子来说明它是如何工作的。但是您似乎在寻求一种无需重新运行模型即可获得预测区间的方法（您的模型是否非常大且计算量大？）。这也可以做到。如何预测——困难的方式（无需重新运行模型） 对于 MCMC 的每次迭代，根据该迭代对协变量值的后验绘制来模拟所需 x 值的响应。所以假设你想预测 X=10 的值。然后如果迭代 1（老化后）的斜率 = 2，截距 = 1，标准差 = 0.5，则从

中绘制 Y 值

Y=rnorm(1, 1+2*10, 0.5)  

然后重复迭代 2、3、4、5... 这些将是您在 X=10 时响应的后验图。 注意：如果您没有监控您的 JAGS 模型中的标准差，那么您很不走运，需要重新拟合模型。

如何预测 - 简单的方法 - 与工作示例 基本思想是插入（到您的数据中）您想要预测其响应的 x 值，以及相关的 y 值 NA。例如，如果您想要 X=10 的预测区间，您只需在数据中包含点 (10, NA)，并为 y 值设置跟踪监视器。

我将 R 中的 JAGS 与 rjags 包一起使用。下面是一个完整的工作示例，首先模拟数据，然后向数据添加一些额外的 x 值，通过 rjags 在 JAGS 中指定和运行线性模型，并总结结果。 Y[101:105] 包含来自 X[101:105] 的后验预测区间的绘制。请注意，Y[1:100] 仅包含 X[1:100] 的 y 值。这些是我们提供给模型的观察数据，它们永远不会随着模型的更新而改变。

library(rjags)
# Simulate data (100 observations)
my.data <- as.data.frame(matrix(data=NA, nrow=100, ncol=2))
names(my.data) <- c("X", "Y")
# the linear model will predict Y based on the covariate X

my.data$X <- runif(100) # values for the covariate
int <- 2     # specify the true intercept
slope <- 1   # specify the true slope
sigma <- .5   # specify the true residual standard deviation
my.data$Y <- rnorm(100, slope*my.data$X+int, sigma)  # Simulate the data

#### Extra data for prediction of unknown Y-values from known X-values
y.predict <- as.data.frame(matrix(data=NA, nrow=5, ncol=2))
names(y.predict) <- c("X", "Y")
y.predict$X <- c(-1, 0, 1.3, 2, 7)

mydata <- rbind(my.data, y.predict)


set.seed(333)
setwd(INSERT YOUR WORKING DIRECTORY HERE)
sink("mymodel.txt")
cat("model

    # Priors

    int ~ dnorm(0, .001)
    slope ~ dnorm(0, .001)
    tau <- 1/(sigma * sigma)
    sigma ~ dunif(0,10) 

    # Model structure

    for(i in 1:R)
    Y[i] ~ dnorm(m[i],tau)
    m[i] <- int + slope * X[i]
    
    ", fill=TRUE)
sink()
jags.data <- list(R=dim(mydata)[1], X=mydata$X, Y=mydata$Y)

inits <- function()list(int=rnorm(1, 0, 5), slope=rnorm(1,0,5),
                         sigma=runif(1,0,10))

params <- c("Y", "int", "slope", "sigma")

nc <- 3
n.adapt <-1000
n.burn <- 1000
n.iter <- 10000
thin <- 10
my.model <- jags.model('mymodel.txt', data = jags.data, inits=inits, n.chains=nc, n.adapt=n.adapt)
update(my.model, n.burn)
my.model_samples <- coda.samples(my.model,params,n.iter=n.iter, thin=thin)
summary(my.model_samples)

【讨论】：

如果您有 50 万个观测值要预测怎么办？有没有办法在不花一年时间的情况下做到这一点“简单的方法”？

你最好还是用“艰难的方式”来做，这并不难。如果您使用 R，您可以对超过 50 万个点的预测进行矢量化处理，这将非常快，并将其包装在一个循环中，该循环遍历 MCMC 链的后迭代。