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

Question

我有以下潜在变量模型:Person j有两个潜在变量，X_j1和X_j2。我们唯一能观察到的是它们的最大值，Y_j = max(X_j1， X_j2)。潜变量为二元正态;它们都有均值mu，方差sigma²，它们的相关性是。我想估计三个参数(mu, sigma²， rho)，只使用Y_j，数据来自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中实现这样的模型(在这种情况下使用了一点hack，因为您正在建模最大值而不是两个变量的和)。但是下面的代码可以编译和运行(尽管它在任何实际意义上都不"工作"——见后面):

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)

有两点需要注意:

1. JAGS不够聪明，无法在满足dsum(max(X[i，]))约束的情况下初始化X的矩阵——所以我们必须使用合理的值来初始化JAGS的X。在这种情况下，我使用的是模拟值，这是作弊-你得到的答案高度依赖于X初始值的选择，而在现实世界中，你不会有模拟值可以依靠。

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

就我所能看到的"我如何编写这个"问题的其他潜在解决方案将落入本质上相同的非混合陷阱，我认为这是这里的一个基本问题(尽管我可能是错的，并且非常欢迎工作的bug/JAGS/Stan代码说明其他情况)。

解决混合失败的方法更难，尽管类似于Carlin &用于模型选择的Chibb方法可能有效(强制min(pseudo_X)参数等于Y以鼓励切换)。这可能是棘手的工作，但如果你能得到一个人的帮助，与bug/JAGS的合理数量的经验，你可以尝试它-见:卡林，b.p.，中国，美国，1995。基于马尔可夫链蒙特卡罗方法的贝叶斯模型选择。J. R.州社爵士。B 57, 473-484 .

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

顺便说一下，我意识到这并不一定回答你的问题，但我认为这是一个有用的例子，说明了"它是可编码的"one_answers"它是可行的"之间的区别。

马特

p。一个更聪明的解决方案是直接考虑两个正态变量的最大值的分布——我不确定这样的分布是否存在，但如果存在，并且你可以得到它的PDF，那么分布可以直接使用0/1技巧进行编码，而不必考虑最小值的值。

Answer 2

我相信您可以在Stan语言中对此进行建模，将可能性视为具有相等权重的两分量混合物。Stan代码看起来像

data {
  int<lower=1> N;
  vector[N] Y;
}
parameters {
  vector<upper=0>[2] diff[N];
  real mu;
  real<lower=0> sigma;
  real<lower=-1,upper=1> rho;
}
model {
  vector[2] case_1[N];
  vector[2] case_2[N];
  vector[2] mu_vec;
  matrix[2,2] Sigma;
  for (n in 1:N) {
    case_1[n][1] = Y[n]; case_1[n][2] = Y[n] + diff[n][1];
    case_2[n][2] = Y[n]; case_2[n][1] = Y[n] + diff[n][2];
  }
  mu_vec[1] = mu; mu_vec[2] = mu;
  Sigma[1,1] = square(sigma);
  Sigma[2,2] = Sigma[1,1];
  Sigma[1,2] = Sigma[1,1] * rho;
  Sigma[2,1] = Sigma[1,2];
  // log-likelihood
  target += log_mix(0.5, multi_normal_lpdf(case_1 | mu_vec, Sigma), 
                         multi_normal_lpdf(case_2 | mu_vec, Sigma));
  // insert priors on mu, sigma, and rho
}