Consequences of putting multiple models in model block

[Quinn Lathrop]

Hi,

In working with simulations, I often have competing models fit to the same data. Each model is kept in a different .stan file, but they often contain very similar syntax. My question here is, what are the consequences of putting multiple models within the same model block? More simply, is there any situation where this is a reasonable approach?

I have found this helpful in communicating the differences in models and in understanding how to implement those models in stan. So the question is, is this a valid use of stan? Or should I never write model code like this ever again?

Example: Mean and Covariance Under Missing Data

Let's say we want to estimate the mean and covariance from two variables x and y. In our sample, we always observe x, but sometimes y is missing. Two basic modeling choices are to use only Complete Cases (only use data where x and y are observed) or to not delete any data and use Available Data methods. Outside of stan, the complete case could be implemented with na.omit(), or cov(use = ""complete.obs") and the available data method through, for example, an EM algorithm as in the package rsem. Of course deleting data is not a good thing and everyone should know that. While this simple example demonstrates this, being in stan, it lends itself to then be expanded into more complicated models (like selection models).

 Here we go:

#In R
library(MASS)
library(rstan)

set.seed(123)

n <- 500
xy_true <- mvrnorm(n, c( 0, 0 ), matrix( c( 1, .8, .8, 1), 2) )
colnames(xy_true) <- c("x", "y")

r <- ifelse(xy_true[,1] > .7, 0, 1)
    #missing indicator, r = 1 obs, r = 0 for missing
    #this is a simple and harsh MAR mechanism (missing on y is determined by x)
    #about 25% missingness on y

xy_obs <- xy_true[ r == 1, ]
n_obs <- nrow(xy_obs)

xy_mis <- xy_true[ r == 0, ]
xy_mis[,2] <- NA #missing on y
x_mis <- xy_mis[,1]
n_mis <- nrow(xy_mis)

sdata <- list("n_obs" = n_obs,
                  "n_mis" = n_mis,
                  "xy_obs" = xy_obs,
                  "x_mis" = x_mis)

sout <- stan(model_code = mod_code, data = sdata, pars = c("mu_1", "sigma_1", "mu_2", "sigma_2"), iter = 400)

Where the model is (notice the two competing models).

mod_code <- 
'

data {

  int<lower=0> n_obs; //number of obs cases

  int<lower=0> n_mis; //number of missing cases

  vector[2] xy_obs[n_obs]; //xy for observed cases

  real x_mis[n_mis]; //x for cases with missing y

}

parameters {
  // Complete Case Parameters

  vector[2] mu_1; 

  cov_matrix[2] sigma_1;

  // Available Data Parameters

  vector[2] mu_2;

  cov_matrix[2] sigma_2;

  real y_mis[n_mis];

}

transformed parameters{

  vector[2] xy_mis[n_mis]; //xy for cases with missing y

  
                     // I could not figure out how to vectorize this loop
                     // seems simple. Any tips?

  for(i in 1:n_mis){

    xy_mis[i,1] <- x_mis[i]; //x_mis is obsvered

    xy_mis[i,2] <- y_mis[i]; //y_mis is parameter

  }

}

// This model block has pedagogical qualities

model {

    //Complete Case Model
    xy_obs ~ multi_normal(mu_1, sigma_1);

    //Available Data Model

    xy_obs ~ multi_normal(mu_2, sigma_2);

    xy_mis ~ multi_normal(mu_2, sigma_2);

}

'

And then print(sout) or launch_shinystan(sout) clearly shows the bias in the mean and covariance under the Complete Case model. To understand why there are differences in the parameter estimates, the answer is in the stan model block. And so forth.

Thanks in advance for any comments,

[Bob Carpenter]

Nothing wrong in principle, but the algorithms tend to be
faster done separately because of coupling all the systems
step-by-step in HMC.

One thing you can try to do is pull out functionality into
functions, but don't (yet) have a way to set up a library of
functions to include (but it's on my to-do list for this year).

More on code inline.

- Bob

> On Mar 16, 2016, at 6:20 PM, Quinn Lathrop <quinn....@gmail.com> wrote:
>
> Hi,
>
> In working with simulations, I often have competing models fit to the same data. Each model is kept in a different .stan file, but they often contain very similar syntax. My question here is, what are the consequences of putting multiple models within the same model block? More simply, is there any situation where this is a reasonable approach?
>
> I have found this helpful in communicating the differences in models and in understanding how to implement those models in stan. So the question is, is this a valid use of stan? Or should I never write model code like this ever again?
>

> ...

> data {
> int<lower=0> n_obs; //number of obs cases
> int<lower=0> n_mis; //number of missing cases
> vector[2] xy_obs[n_obs]; //xy for observed cases
> real x_mis[n_mis]; //x for cases with missing y
> }
>
> parameters {
> // Complete Case Parameters
> vector[2] mu_1;
> cov_matrix[2] sigma_1;
>
> // Available Data Parameters
> vector[2] mu_2;
> cov_matrix[2] sigma_2;
> real y_mis[n_mis];
> }

You really want to parameterize this in terms of Cholesky
factors. We strongly recommend scaling a Cholesky factor for a
correlation matrix (with an LKJ prior), but if you really want
the improper uniform on covariance matrices, it's
much more arithmetically stable and efficient to use the Cholesky
forms (of both variables and densities).

> transformed parameters{
> vector[2] xy_mis[n_mis]; //xy for cases with missing y
>
> // I could not figure out how to vectorize this loop
> // seems simple. Any tips?
> for(i in 1:n_mis){
> xy_mis[i,1] <- x_mis[i]; //x_mis is obsvered
> xy_mis[i,2] <- y_mis[i]; //y_mis is parameter
> }
> }

You can do this:

x[, 1] <- x_mis;
x[, 2] <- y_mis;

> // This model block has pedagogical qualities

I'd avoid this kind of comment, too --- comments overall make
code *less* readable if the code is reasonable.

> model {
>
> //Complete Case Model
> xy_obs ~ multi_normal(mu_1, sigma_1);
>
> //Available Data Model
> xy_obs ~ multi_normal(mu_2, sigma_2);
> xy_mis ~ multi_normal(mu_2, sigma_2);

Even more efficient if you append_put xy_mis together into
a single array of vectors. But in this case, I don't think
the xy_mis ~ ... line will affect the estimates of mu_2 and
sigma_2; instead, you could probably generate the predictions
for missing y in generated quantities more efficiently (with
help from someone better than me at conditioning a multivariate
normal). But what you have will also work.

- Bob