r - 当观察到的节点是最大潜在节点时使用 JAGS 或 STAN

Question

我有以下潜变量模型:人 j 有两个潜变量，Xj1 和 Xj2。我们唯一能观察到的是它们的最大值，Yj = max(Xj1, Xj2)。潜在变量是二元正态的；它们每个都有均值 mu、方差 sigma2，并且它们的相关性是 rho。我想仅使用 Yj 来估计三个参数(mu、sigma2、rho)，数据来自 n 个患者，j = 1,...,n。

我试图在 JAGS 中拟合这个模型(所以我将先验放在参数上)，但我无法编译代码。这是我用来调用 JAGS 的 R 代码。首先，给定参数的一些真实值，我生成数据(潜在变量和观察变量):

# true parameter values
mu <- 3
sigma2 <- 2
rho <- 0.7

# generate data
n <- 100
Sigma <- sigma2 * matrix(c(1, rho, rho, 1), ncol=2)
X <- MASS::mvrnorm(n, c(mu,mu), Sigma) # n-by-2 matrix
Y <- apply(X, 1, max)

然后我定义了 JAGS 模型，并编写了一个小函数来运行 JAGS 采样器并返回样本:

# JAGS model code
model.text <- '
model {
  for (i in 1:n) {
    Y[i] <- max(X[i,1], X[i,2]) # Ack!
    X[i,1:2] ~ dmnorm(X_mean, X_prec)
  }

  # mean vector and precision matrix for X[i,1:2]
  X_mean <- c(mu, mu)
  X_prec[1,1] <- 1 / (sigma2*(1-rho^2))
  X_prec[2,1] <- -rho / (sigma2*(1-rho^2))
  X_prec[1,2] <- X_prec[2,1]
  X_prec[2,2] <- X_prec[1,1]

  mu ~ dnorm(0, 1)
  sigma2 <- 1 / tau
  tau ~ dgamma(2, 1)
  rho ~ dbeta(2, 2)
}
'

# run JAGS code. If latent=FALSE, remove the line defining Y[i] from the JAGS model
fit.jags <- function(latent=TRUE, data, n.adapt=1000, n.burnin, n.samp) {
  require(rjags)
  if (!latent)
    model.text <- sub('\n *Y.*?\n', '\n', model.text)
  textCon <- textConnection(model.text)
  fit <- jags.model(textCon, data, n.adapt=n.adapt)
  close(textCon)
  update(fit, n.iter=n.burnin)
  coda.samples(fit, variable.names=c("mu","sigma2","rho"), n.iter=n.samp)[[1]]
}

最后，我调用 JAGS，只提供观察到的数据:

samp1 <- fit.jags(latent=TRUE, data=list(n=n, Y=Y), n.burnin=1000, n.samp=2000)

遗憾的是，这会导致错误消息:“Y[1] 是一个逻辑节点，无法观察”。 JAGS 不喜欢我使用“<-”为 Y[i] 赋值(我用“Ack!”表示违规行)。我理解投诉，但我不确定如何重写模型代码来解决这个问题。

此外，为了证明其他一切(除了“Ack!”行)都很好，我再次运行模型，但这次我向它提供 X 数据，假装它确实被观察到了。这运行完美，我得到了很好的参数估计:

samp2 <- fit.jags(latent=FALSE, data=list(n=n, X=X), n.burnin=1000, n.samp=2000)
colMeans(samp2)

如果您能找到一种在 STAN 而不是 JAGS 中对该模型进行编程的方法，那对我来说没问题。

Answer 1

从理论上讲，您可以使用 dsum 分布在 JAGS 中实现这样的模型(在这种情况下，在建模两个变量的最大值而不是总和时使用了一些技巧)。但是以下代码确实可以编译并运行(尽管它在任何真正意义上都不起作用 - 见下文):

set.seed(2017-02-08)

# true parameter values
mu <- 3
sigma2 <- 2
rho <- 0.7

# generate data
n <- 100
Sigma <- sigma2 * matrix(c(1, rho, rho, 1), ncol=2)
X <- MASS::mvrnorm(n, c(mu,mu), Sigma) # n-by-2 matrix
Y <- apply(X, 1, max)

model.text <- '
model {
  for (i in 1:n) {
    Y[i] ~ dsum(max_X[i])
    max_X[i] <- max(X[i,1], X[i,2])
    X[i,1:2] ~ dmnorm(X_mean, X_prec)
    ranks[i,1:2] <- rank(X[i,1:2])
    chosen[i] <- ranks[i,2]
  }

  # mean vector and precision matrix for X[i,1:2]
  X_mean <- c(mu, mu)
  X_prec[1,1] <- 1 / (sigma2*(1-rho^2))
  X_prec[2,1] <- -rho / (sigma2*(1-rho^2))
  X_prec[1,2] <- X_prec[2,1]
  X_prec[2,2] <- X_prec[1,1]

  mu ~ dnorm(0, 1)
  sigma2 <- 1 / tau
  tau ~ dgamma(2, 1)
  rho ~ dbeta(2, 2)

  #data# n, Y
  #monitor# mu, sigma2, rho, tau, chosen[1:10]
  #inits# X
}
'

library('runjags')

results <- run.jags(model.text)
results
plot(results)

有两点需要注意:

JAGS 不够聪明，无法在满足 dsum(max(X[i,])) 约束的同时初始化 X 的矩阵 - 因此我们必须使用合理的值为 JAGS 初始化 X。在这种情况下，我使用的是作弊的模拟值 - 您得到的答案高度依赖于 X 初始值的选择，而在现实世界中，您将无法依靠模拟值。

max() 约束会导致我无法在一般框架内找到解决方案的问题:与通常的 dsum 约束不同，dsum 约束允许一个参数减少而另一个参数增加，因此这两个参数始终被使用，min( ) X[i,] 的值被忽略，因此采样器可以随心所欲地做。这将很少(即永远不会)导致 min(X[i,]) 的值恰好与 Y[i] 相同，这是采样器在两个 X[i] 之间“切换”所需的条件,]。因此切换永远不会发生，并且在初始化时选择作为最大值的 X[] 保持为最大值 - 我添加了一个跟踪参数“选择”来说明这一点。

据我所知，“我如何编码这个”问题的其他潜在解决方案将陷入本质上相同的非混合陷阱，我认为这是这里的一个基本问题(尽管我可能是错的，并且非常欢迎工作BUGS/JAGS/Stan 代码，否则说明)。

混合失败的解决方案更难，尽管类似于用于模型选择的 Carlin & Chibb 方法可能有效(强制 min(pseudo_X) 参数等于 Y 以鼓励切换)。这可能很难开始工作，但是如果您可以从具有合理 BUGS/JAGS 经验的人那里获得帮助，您可以尝试一下 - 请参阅:
Carlin, B.P., Chib, S., 1995. 通过马尔可夫链蒙特卡罗方法选择贝叶斯模型。 J. R. 统计社会。爵士。 B 57, 473–484。

或者，您可以尝试稍微不同地考虑这个问题，并将 X 直接建模为一个矩阵，其中第一列全部缺失，第二列全部等于 Y。然后您可以使用 dinterval() 对它们的缺失值设置约束必须低于相应的最大值。我不确定这在估计 mu/sigma2/rho 方面效果如何，但可能值得一试。

顺便说一句，我意识到这不一定能回答您的问题，但我认为这是“是否可编码”和“是否可行”之间区别的一个有用示例。

马特

附:一个更聪明的解决方案是直接考虑两个正态变量的最大值的分布 - 我不确定这样的分布是否存在，但它确实存在，您可以获得它的 PDF，然后可以直接使用分布对分布进行编码零/一技巧而无需考虑最小值的值。