2

我想使用 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)
4

1 回答 1

4

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

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

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)
于 2015-11-15T17:03:26.200 回答