Chains getting stuck in a hierarchical nested domain model

[Matt]

Hi Stan!
I've been implementing model similar to the one described by Thurston et al. 2009 (Bayesian models for multiple outcomes nested in domains - http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3031784/). It shrinks individual outcomes towards a common domain effect, and each domain towards a common overall effect. I had a good JAGS implementation until I ran into dependent variables that correlated quite highly (r > ~ 0.8) and this severely decreased the number of effective samples. So I implemented a similar model in STAN (minimal example below), but I'm still getting a low number of effective samples which seem to be due to the chains getting sporadically stuck.

I've tried re-parameterising with the Matt trick, using "nondiag_mass=T", and adjusting/increasing delta. Sometimes I get a really good number of effective samples, albeit after a longer wait, but then when I make slight adjustments to delta (up or down), it can go back to being stuck. This means I have to keep fiddling to get nice settings.

Is there a better way to re-parameterise this model to improve chain stickiness without compromising speed?
Short of that, it's easy enough to simply remove the variable as it's not supplying much information, and that's what I'll prob do atm.

Thanks,
Matt

require(rstan)
Ntot  <- 50                   #Sample size
dom   <- c(1,1,1,2,2)         #Outcome domains indicators
Nvar  <- length(dom)
x     <- rnorm(Ntot, 0, 1)
y     <- sapply(1:Nvar, function(i) x + rnorm(Ntot, 0, dom[i]))
y[,3] <- y[,1] + y[,2]        #Create highly correlated variables

zx       <- scale(x)
zy       <- apply(y, 2, scale)
Ndom     <- length(unique(dom))
dataList <- list(x=zx[,1], y=t(zy), Nvar=Nvar, Ntot=Ntot, id=1:Ntot, dom=dom, Ndom = Ndom)

mstr<- '
data {
  int Ntot;
  int Nvar;
  int Ndom;
  int dom[Nvar];
  int id[Ntot];
  vector[Ntot] x;
  matrix[Nvar, Ntot] y;
}
parameters {
  real<lower=0,upper=10> sd[Nvar];  //SD per y variable
  real b1;                          //Overarching beta
  vector[Nvar] e1b;                 //MTrick for Ind Outcome betas
  real<lower=0> b1sd;
  vector[Ndom] e1d;                 //MTrick for Domain betas
  real<lower=0> b1dsd;
  vector[Ntot] e0bs;                //MTrick for participant random intercepts
  real<lower=0,upper=10> sd0s;      
  matrix[Ndom, Ntot] e0bds;         //MTrick for participant:Domain random intercepts
  real<lower=0,upper=10> sd0d[Ndom];
}
transformed parameters{
  vector[Nvar] B1O;                 //Output for each individual beta
  vector[Nvar] b1o;                 //Outcome betas
  vector[Ndom] b1d;                 //Domain betas
  vector[Ntot] b0s;                 //Participant intercepts
  matrix[Ndom, Ntot] b0ds;          //Participant:Domain intercepts
  matrix[Nvar, Ntot] mu;

  b1o <- e1b * b1sd;
  b1d <- e1d * b1dsd;

  b0s  <- e0bs * sd0s;
  for(d in 1:Ndom){
    b0ds[d] <- e0bds[d] * sd0d[d];
  }

  for(np in 1:Nvar){
    for(i in 1:Ntot){
      mu[np, i] <- b0s[id[i]] + b0ds[dom[np], id[i]] + (b1 + b1d[dom[np]] + b1o[np]) * x[i];
    }
  }
//Output
  for(np in 1:Nvar){
    B1O[np] <- b1 + b1d[dom[np]] + b1o[np];
  }
}
model {
  b1    ~ normal(0, 10);
  e1b   ~ normal(0, 1);
  e1d   ~ normal(0, 1);
  e0bs  ~ normal(0, 1);
  for(d in 1:Ndom){
    e0bds[d] ~ normal(0, 1);
  }
  b1sd  ~ cauchy(0, 10);
  b1dsd  ~ cauchy(0, 10);
// Likelihood
  for(np in 1:Nvar){
    y[np] ~ normal(mu[np], sd[np]);
  }
}
'
seed = 69
f1 <- stan(model_code = mstr, data = dataList, chains = 4, iter=100,
           seed = seed, chain_id = 1,# nondiag_mass = T,# delta = 0.9,
           pars = c("b1", "b1o", "b1d", "sd", "b1sd", "b1dsd", "sd0d", "B1O"))

f2 <- stan(fit=f1, data = dataList, chains = 4, iter=2500,
           seed = seed, chain_id = 1, #nondiag_mass = T, delta = 0.9,
           pars = c("b1", "b1o", "b1d", "sd", "b1sd", "b1dsd", "sd0d", "B1O"))

[Bob Carpenter]

Although it won't help with n_eff, for speed, you definitely
want to vectorize.

for(d in 1:Ndom)
e0bds[d] ~ normal(0, 1);

should be written as:

e0bds ~ normal(0,1);

Unfortunately, this

for(np in 1:Nvar)
y[np] ~ normal(mu[np], sd[np]);

is as far as you can vectorize because mu is a
matrix. It'd be a bit more efficient if mu and y were
arrays of vectors because it eliminates copies.

You should move what you called "output"
out of transformed parameters to a generated quantities
block after the model:

generated quantities {

vector[Nvar] B1O; //Output for each individual beta

for(np in 1:Nvar)
B1O[np] <- b1 + b1d[dom[np]] + b1o[np];
}

This will save a lot of unnecessary computation.


It won't help as much, but you can also replace some of
the loops with their matrix equivalents.

for(d in 1:Ndom)
b0ds[d] <- e0bds[d] * sd0d[d];

can use elementwise operations:

b0ds <- e0bds .* sd0d;

If you're having trouble with parameters wandering too
far afield, you can tighten the priors, too.

You should also put appropriate upper and lower-bound
constraints on the data variables dom and id for clarity
and data checking --- it won't change how things run if
everything's valid.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Matt]

Thanks Bob,
I'll try replacing all matrix objects to vector objects (I found the matrix notations easier to track initially) and move the output to the end.

Just to clarify, "e0bds" is a matrix[Ndom, Ntot]. Will "e0bds ~ normal(0,1)" sample for the whole matrix if I replace it with an equivalent vector?

Matt

[Bob Carpenter]

On 7/3/13 1:39 PM, Matt wrote:
> Thanks Bob,
> I'll try replacing all matrix objects to vector objects (I found the matrix notations easier to track initially) and
> move the output to the end.

With matrix objects, it's a bit more efficient to pull
out by column (even though the indexing notation only
supports rows). That's because the underlying matrix lib
we use, Eigen, organizes matrices by column (by default,
which is what we went with). Arrays of row vectors
(or arrays of vectors) are more efficient because there's
no copying needed.

> Just to clarify, "e0bds" is a matrix[Ndom, Ntot]. Will "e0bds ~ normal(0,1)" sample for the whole matrix if I replace it
> with an equivalent vector?

No. Sorry. I didn't realize that one was a matrix.

We'll get there eventually, but it doesn't work yet.
But it's not that much overhead to call once per row,
especially if the rows are long.

- Bob