1

我是 JAGS 的新手，我正在尝试使用 9 个非连续预测变量来预测二元结果（0/1）。预测值可能是 0、1 或 2。这是我第一次这样做，尽管我可以让模型运行，但我 100% 肯定这里肯定存在许多问题。

数据文件样本（列表）

$y
[1] 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0

$N
[1] 50

$oAnt
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1
[29] 1 1 1 1 1 1 2 1 1 0 1 1 1 1 1 2 1 1 1 1 1 1

$nAnt
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

$cAnt
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0
[29] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0

$oPen
[1] 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1 1 1 0
[29] 1 0 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 0 2 1 1

$nPen
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1 2 1
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1

$cPen
[1] 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0
[29] 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0

$oFin
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1
[29] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

$nFin
[1] 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1 1 1 2 1 1 1
[29] 1 1 1 1 2 1 1 1 1 2 1 1 1 2 1 1 1 1 1 1 1 3

$cFin
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1
[29] 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0

模型

model {
    for( i in 1 : N ){
        y[i] ~ dbern(mu[i])
        mu[i] <- 1/(1+exp(-(b0 + b1*oAnt[i] + b2*nAnt[i] + b3*cAnt[i] + b4*oPen[i] + b5*nPen[i] + b6*cPen[i] + b7*oFin[i] + b8*nFin[i] + b9*cFin[i])))
}
b0 ~ dnorm(0, 1.0e-12)
b1 ~ dnorm(0, 1.0e-12)
b2 ~ dnorm(0, 1.0e-12)
b3 ~ dnorm(0, 1.0e-12)
b4 ~ dnorm(0, 1.0e-12)
b5 ~ dnorm(0, 1.0e-12)
b6 ~ dnorm(0, 1.0e-12)
b7 ~ dnorm(0, 1.0e-12)
b8 ~ dnorm(0, 1.0e-12)
b9 ~ dnorm(0, 1.0e-12)
}

我使用来自glm()模型的估计值作为初始值（正如 A. Gelman 所建议的那样）——但为了简单起见，我们假设我会让 JAGS 选择链的初始值。

跑步模型等

jagsModel = jags.model(file = "antPenFin.txt", data = dataList, n.chains = 2, n.adapt = 500)

update(jagsModel, n.iter = 500)

codaSamples = coda.samples(jagsModel,
variable.names = names(dataList)[3:11], n.iter = 5000)

问题

我的模型的输出看起来完全不正确（当我尝试绘制它时会变得很清楚）。我敢肯定这里有一些非常基本的问题。有人可以帮忙吗？

非常感谢。

4

1 回答 1

0

您的问题是您无法提供一系列数字（在您的情况下为 0 到 3）作为分类协变量。目前，您的模型将这些数字解释为连续的。您需要做的是将这些转换为虚拟变量。您可以使用model.matrixR 中的函数轻松完成此操作。我将生成一些数据作为示例，然后您可以将其应用于您的数据。

# generate y
y <- sample(0:1, 30, replace = TRUE)
# add three different categorical covariates. Note here that these are all
# factors.
oAnt <- factor(sample(0:2, 30, replace = TRUE))
cAnt <- factor(sample(0:1, 30, replace = TRUE))
nFin <- factor(sample(0:3, 30, replace = TRUE))

# create your model matrix
my_matrix <- model.matrix(y ~ oAnt + cAnt + nFin)

head(my_matrix)

   (Intercept) oAnt1 oAnt2 cAnt1 nFin1 nFin2 nFin3
1           1     1     0     1     0     1     0
2           1     1     0     1     1     0     0
3           1     1     0     0     0     0     1
4           1     1     0     1     0     0     1
5           1     0     0     1     0     1     0
6           1     1     0     1     0     0     1

对于预测变量的因子中的每个级别，您将需要创建 n-1 列（例如，nFin 范围为 0-3，您将创建 3 列虚拟变量）。因此，您将在模型中包含更多回归系数。制作模型矩阵后，您可以将其转换为列表。

# remove the intercept from the list as you already have it in your model
matrix_list <- as.list(data.frame(my_matrix[,-1]))

从那里，您需要做的就是为模型创建更多预测变量。此外，如果nAnt您的模型中确实是所有的，请继续删除它，因为您基本上只包括两个截距。

于 2016-06-16T13:46:44.433 回答