Using JAGS or STAN when an observed node is the max of latent nodes

Question

I have the following latent variable model: Person j has two latent variables, X_j1 and X_j2. The only thing we get to observe is their maximum, Y_j = max(X_j1, X_j2). The latent variables are bivariate normal; they each have mean mu, variance sigma², and their correlation is rho. I want to estimate the three parameters (mu, sigma², rho) using only Y_j, with data from n patients, j = 1,...,n.

I've tried to fit this model in JAGS (so I'm putting priors on the parameters), but I can't get the code to compile. Here's the R code I'm using to call JAGS. First I generate the data (both latent and observed variables), given some true values of the parameters:

# 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)

Then I define the JAGS model, and write a little function to run the JAGS sampler and return the samples:

# 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]]
}

Finally, I call JAGS, feeding it only the observed data:

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

Sadly this results in an error message: "Y[1] is a logical node and cannot be observed". JAGS does not like me using "<-" to assign a value to Y[i] (I denote the offending line with an "Ack!"). I understand the complaint, but I'm not sure how to rewrite the model code to fix this.

Also, to demonstrate that everything else (besides the "Ack!" line) is fine, I run the model again, but this time I feed it the X data, pretending that it's actually observed. This runs perfectly and I get good estimates of the parameters:

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

If you can find a way to program this model in STAN instead of JAGS, that would be fine with me.

Answer 1

Theoretically you can implement a model like this in JAGS using the dsum distribution (which in this case uses a bit of a hack as you are modelling the maximum and not the sum of the two variables). But the following code does compile and run (although it does not 'work' in any real sense - see later):

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)

Two things to note:

1. JAGS isn't smart enough to initialise the matrix of X while satisfying the dsum(max(X[i,])) constraint on its own - so we have to initialise X for JAGS using sensible values. In this case I'm using the simulated values which is cheating - the answer you get is highly dependent on the choice of initial values for X, and in the real world you won't have the simulated values to fall back on.

2. The max() constraint causes problems to which I can't think of a solution within a general framework: unlike the usual dsum constraint that allows one parameter to decrease while the other increases and therefore both parameters are used at all times, the min() value of X[i,] is ignored and the sampler is therefore free to do as it pleases. This will very very rarely (i.e. never) lead to values of min(X[i,]) that happen to be identical to Y[i], which is the condition required for the sampler to 'switch' between the two X[i,]. So switching never happens, and the X[] that were chosen at initialisation to be the maxima stay as the maxima - I have added a trace parameter 'chosen' which illustrates this.

As far as I can see the other potential solutions to the 'how do I code this' question will fall into essentially the same non-mixing trap which I think is a fundamental problem here (although I might be wrong and would very much welcome working BUGS/JAGS/Stan code that illustrates otherwise).

Solutions to the failure to mix are harder, although something akin to the Carlin & Chibb method for model selection may work (force a min(pseudo_X) parameter to be equal to Y to encourage switching). This is likely to be tricky to get working, but if you can get help from someone with a reasonable amount of experience with BUGS/JAGS you could try it - see: Carlin, B.P., Chib, S., 1995. Bayesian model choice via Markov chain Monte Carlo methods. J. R. Stat. Soc. Ser. B 57, 473–484.

Alternatively, you could try thinking about the problem slightly differently and model X directly as a matrix with the first column all missing and the second column all equal to Y. You could then use dinterval() to set a constraint on the missing values that they must be lower than the corresponding maximum. I'm not sure how well this would work in terms of estimating mu/sigma2/rho but it might be worth a try.

By the way, I realise that this doesn't necessarily answer your question but I think it is a useful example of the difference between 'is it codeable' and 'is it workable'.

Matt

ps. A much smarter solution would be to consider the distribution of the maximum of two normal variates directly - I am not sure if such a distribution exists, but it it does and you can get a PDF for it then the distribution could be coded directly using the zeros/ones trick without having to consider the value of the minimum at all.

Answer 2

I believe you can model this in the Stan language treating the likelihood as a two component mixture with equal weights. The Stan code could look like

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
}