Sudden drops in "accpet_stat__" and "treedepth__"

[Feezdev]

Hi everyone,

I have a Stan model which is performing a hierarchical logistic regression on about 100 subjects, each performing multiple trials.

I think I might be having some convergence issues as it seems that the model is sampling fine on up to a point where the "accept_stat__" drops suddenly and the chains then start to make very small moves (see attached traceplots for hyper-parameters and diagnostics -- gray area is warmup time).

I'm trying to understand why this is going on by trying to identify if some parameter gets stuck away from a "typical value" when the model is sampling from in the areas where the "accept_stat__" and "treedepth__" diagnostics drop (but can't seem to find anything out of the ordinary) :-/

I've attached the traceplots of 2 different runs (each with 2 chains) on the exact same model and same data.

Here is the model (if it helps):

data{

    int<lower=1> Ns;

    int<lower=0> Nx;

    int<lower=1> Ni;

    int<lower=0,upper=1> y[Ns,Nx];

    matrix<lower=0>[Ns,Nx] x1;

    matrix<lower=0>[Ns,Nx] x2;

}

parameters{

   //Grp Lvl

   real<lower=0> mu_intercept;

   real<lower=0> a_guess;

   real<lower=0> mu_x1;         // expected to be positive

   real mu_x2;                  // expected to be negative

   real mu_x2quad;              // expected to be negative

   real<lower=0.01> b_guess;

   real<lower=0.01> sd_intercept;

   real<lower=0.01> sd_x1;

   real<lower=0.01> sd_x2;

   real<lower=0.01> sd_x2quad;

   //Subj Lvl

   vector<lower=0, upper=20>[Ns] theta_intercept;

   vector<lower=0, upper=1 >[Ns] theta_guess;

   vector<lower=0,upper=15 >[Ns] theta_x1;

   vector<lower=-50,upper=1>[Ns] theta_x2; 

   vector<lower=-50,upper=1>[Ns] theta_x2quad;

}

model{

  mu_intercept ~ normal(6,100);

  mu_x1 ~ normal(6,100);

  mu_x2 ~ normal(-25,100);

  mu_x2quad ~ normal(-25,100);

  sd_intercept ~ cauchy(0,10);

  sd_x1 ~ cauchy(0,10);

  sd_x2 ~ cauchy(0,25);

  sd_x2quad ~ cauchy(0,25);

  theta_intercept ~ normal(mu_intercept,sd_intercept);

  theta_x1 ~ normal(mu_x1,sd_x1);

  theta_x2 ~ normal(mu_x2,sd_x2);

  theta_x2quad ~ normal(mu_x2quad,sd_x2quad);

  a_guess ~ normal(1,1);

  b_guess ~ normal(9,3);

  theta_guess ~ beta(a_guess,b_guess);

  for (i_subj in 1:Ns) {

    for (i_trial in 1:Nx) {

      y[i_subj,i_trial] ~ bernoulli( (theta_guess[i_subj]*(0.5))

                                  +  ( (1-theta_guess[i_subj]) * inv_logit( theta_intercept[i_subj]

                                                                     + x1[i_subj,i_trial]*theta_x1[i_subj]

                                                                     + (x2[i_subj,i_trial])*theta_x2[i_subj]

                                                                     + (x2[i_subj,i_trial]^2)*theta_x2quad[i_subj]

                                                                     )

                                     )

                                   );

    }

  }

}

I can't seem to understand why the model would sample fine on one run, then stop sampling well half way through on run1. Then on the same model/data on run2 it only happens "sporadically".

Is this bad (my guess is yes)? Am I doing something very wrong here?

Best,
Fd

[Feezdev]

Dear Stan team/users,

I'm updating the previous post as I think I've managed to fix my issue and thought that it could be useful to other beginners. (tldr: used the "non-centered" Matt trick)

My first mistake was to assume that because an earlier/simpler model could successfully recover parameters from surrogate data, then I could continue building on more complicated models without further checks and that I'd be ok (how wrong I was!).

To attempt to fix my problem: 

I regenerated surrogate data for that more complex model, then tried to reparametrize the model using different alternatives until it worked (thankfully it worked before I ran out of alternative reparametrizations):

First I tried to rescale the inputs by normalizing them (z-scored). It seems to have helped a bit as the chains were mixing a bit better than before but still eventually got stuck sometimes with low "accept_stat__" & "tree_depth__". 

Then I attempted to use the non-centered (Matt Trick): 

1) First I tried the Matt trick on subject-level parameters while keeping the initial hyperpriors on different scale/locations. That worked okay-ish, but not perfectly.

2) Then tried the Matt trick but removed nearly all lower and upper constraint on parameters, all parameters centered on zero with unit variance. -- Mixing was now working really well and no divergent traces! However parameters recovered were not as good as these recovered from previous attempt with hyperparameters on different scales/locations (my assumption then was that my hyperpriors were too tight).

3) Tried Matt trick, again without upper lower constraints on parameters, centered on zero but with wider variance (i.e. N(0,10) & Cauchy(0,2.5)). -- This worked like a charm in terms of both mixing & parameter recovery (K-ching!) 

4) Then fit the model to real data, and the mixing issues were gone there too (attached) :)

In retrospect, I learnt that it'd always be worth checking that the most complex model from a set of models can also recover parameters from surrogate data, then simpler models deriving from the complex one should also be ok. 

I guess my mistake was to use the "start simple & build up" motto too literally, and I then forgot to check whether things got wrong along the way.

---

Now however, I'm still left relatively puzzled as to how/why this trick works at all!? (Please forgive me if this is really naive)

I feel like by using Matt's trick, I'm only moving the scale/location problem from the model block up to the transformed parameter block, but the parameters are essentially the same then?!

Also: Why is "non-centered" Matt-trick called "non-centered"? Maybe this is related to my previous question about not understanding what the Matt trick does, as once reparametrized, the parameters are now all at the same location in the model_block (i.e. "centered" on zero)? 

Finally, I've seen other scripts using Matt-trick together with Phi_approx() or Phi() in the transformed parameter block. 

I assume to force a parameter to be in between [0 and 1] or [0 and any Number] by adding a constant (please correct me if I'm wrong). 

But what would be the benefit of, say, using N(0,1) then -> phi_approx( mu + sd * N(0,1)) vs. using directly a Beta distribution in the model block?

Also just wanted to say thank you so much to the community, Stan is amazing!

Best,

[Bob Carpenter]

You want to use a non-centered parameterization.

We also recommend removing all those interval constraints
and just using soft constraints---the intervals cause problems
both statistically (bias) and computationally (Hamiltonians running
to infinity).

- 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.
> <Hyperparams_run1.png><Hyperparams_run2.png><Diagnostic_Run1.png><Diagnostic_Run2.png>

[Michael Betancourt]

https://arxiv.org/abs/1312.0906

https://github.com/betanalpha/knitr_case_studies/tree/master/divergences_and_bias

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

<Diag_non_centered.png>

[Feezdev]

Thanks a lot for the links to the publications Michael, I'll have a look :)

Bob, thanks a lot for the advice! That's exactly what I ended up finding out as well by tinkering with the model and fake data (see previous post) :)