Low n_eff in a hierarchical shrinkage problem

[William Thomson]

Hi all,

I've been working on a problem as part of an industrial placement in which the task is to do some variable selection/shrinkage estimation in a simple nonlinear model. I've coded up a few different shrinkage priors using the scale-mixture framework (Bayesian Lasso, normal-gamma, hierarchical shrinkage/horseshoe) and over the past week or so it's been working quite nicely using an adaptation of the model I'd been given at the beginning. However, the people I'm working for are keen for me to fit a model with some issues. I'm getting really low effective sample sizes and the corresponding high autocorrelation. The problem is basically a nonlinear mixed effects model with covariate effects and looks something like:

p * (1 - exp(-S .* time)) + (1 - p) * (K(S) - exp(-a * S .* time))

with S depending additively on 12 (continuous) covariates with interaction terms. There are 76 parameters in S + 2 structural parameters. They want to include interactions of the covariates with time in S too, this turns the mean response into a mixture of sigmoid functions.

I know that mixtures of exponentials are notoriously tricky but even fixing the constant parameter 'a' causes problems. I'm shrinking the coefficients of S, and have a beta prior on 'p'. K is a saturating function of S with one extra parameter that gets a broad log-normal prior. I'm also taking a differentiable piecewise approximation to the maximum of S and 0 to prevent S becoming negative and causing problems -- I insert a piece of a smooth function near zero and match the slopes either side, is this even worth doing?

I have 169 data points -- the time resolution is not great, I have at most 4 time points for a given set of covariates, and measurements only at 4 distinct times across the whole dataset. It seems this is going to cause an identifiability problem for a model with two exponentials?

I've tried normalising the covariates, observations and time, I'm using (what I think are?) reasonable reparameterisations to get everything as close to the unit scale as I can, and most of the parameters in the parameters block come back as being approximately normal and of O(1) size, with some of the 'penalised' coefficients coming back close to zero.

I'm not seeing divergent transitions for the most part, and if I do it's generally only a couple. Is this a case of the model not being identified or do I have another issue? I've noticed that setting tau_tau very small helps somewhat, but my n_eff is still double digits after 2000 iterations, whereas for simplified versions of the model (e.g. a single exponential or an x/(1+x) type function), I generally get n_eff/n close to 1.

I remember seeing Michael Betancourt saying something about a restriction on the posterior variance depending on the dimensionality of the parameter space, might this be at play? My high-dimensional geometric intuition is not particularly strong I'm afraid!

I'd really like to help them with this and I think they have been working with this model for a while, so if there are any suggestions for ways to get some reasonable inferences out I would really appreciate it! I've pasted my stan code below for a case where I'm using a normal-gamma prior. I'm using the github dev version of rstan on windows 10, with adapt_delta = 0.99, I haven't tinkered with treedepth or any other of the control settings yet. Unfortunately I can't provide the data.

All the best,

Will

functions {

  vector ew_max(vector x, real y){

  vector[num_elements(x)] out;

    

  for (i in 1:num_elements(x)){

  if (x[i]<0)

  out[i]=0;

  else if (x[i]<0.01)

  out[i] = (2*0.01*x[i]^2) / (0.01^2 + x[i]^2);

  else

  out[i] = x[i];

  }

  return out;

  }

}

data{

  int<lower=1> npar;

  int<lower=1> ncov;

  int<lower=1> N;

  vector[N] Time;

  vector[N] yObs;

  matrix[N,npar] modmat;

  real<lower=0> nu;

  real<lower=0> tau_tau;

  matrix[N,ncov] X;

}

parameters{

  vector[npar] theta_raw;

  real<lower=0> lambda;

  real<lower=0> tau;

  real<lower=0,upper = pi()/2> sigma_unif;

  real<lower=0> E_raw;

  real<lower=0, upper=1> fe;

}

transformed parameters{

 vector[npar] theta;

 real<lower=0> sigma;

 real<lower=0> E;

 theta = theta_raw * lambda;

 sigma = tan(sigma_unif);

 E = exp(E_raw+log(1000));

}

model{

 vector[N] S;

 vector[N] J ;

 vector[N] myExponent ;

 vector[N] K ;

 vector[N] Y;

 S = modmat*theta;

 J =  ew_max(S,0);

 myExponent =  J .* X[,11] .* X[,12] .* Time;

 K =  E * J .* X[,11]  ./ (1 + E * J .* X[,11]);

 Y = -fe * expm1(-100 * myExponent)+ (fe-1) * K .* expm1(-myExponent);

 tau ~ exponential(tau_tau);

 lambda ~ gamma(tau,nu);

 theta_raw ~ normal(0,1);

 E_raw ~ normal(0,1);

 fe ~ beta(2,18);

 yObs ~ normal(Y,sigma);

}

[William Thomson]

There's a typo in the parameters block, E_raw shouldn't have a lower bound. But the problem persists unfortunately.

[Bob Carpenter]

I couldn't follow the model details, but all of those
exp() and expm1() calls are very prone to overflow/rounding to
zero.

It's not really a mixture model in the usual sense, just a model
using interpolation to define the mean of a normal (fe is in (0, 1)):

Y = -fe * expm1(-100 * myExponent)+ (fe-1) * K .* expm1(-myExponent);

...
yObs ~ normal(Y,sigma);


Are you hitting the max tree depth? If so, it's probably slow, but
that could explain low n_eff.

You need to dive into the posterior you get and see what shapes you get
out of pairs() --- it scatterplots parameters against each other.
What you want to avoid are banana and funnel shapes and high correlations.

- 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[William Thomson]

Hi Bob,

Thanks for getting back. 

The priors I'm using om the theta [i] are like continuous approximations to a spike and slab, so there are some differentiability issues around zero which seem to have been sorted out by reparameterising in the past, but that trick isn't working now.

Treedepth seems like the big problem. I wasn't getting any warnings with the default setting, but if I manually increase it to 15 almost every iteration hits the max and shinystan goes yellow. I'll try to play around with some settings/brute force it and see if I can improve things, or at least get a better idea of what's going on. Thinning should help, right?

Many thanks!

Will

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/wVbCa7z2J98/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

No, thinning won't help anything other than cut
down on the size of the stored sample and reduce
n_eff a bit. If treedepth is that bad, you probably have
some very pathological places in your prior for the step
size to have to get so small that you need 2^15 leapfrog
steps before making a U-turn. Plotting pairs of
parameters can help isolate the bad geometries (funnel shapes,
banana shapes and very strong correlations).

- Bob

> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/wVbCa7z2J98/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --

> 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.