What to look for in pairs() diagnostic plots?

[David]

I'm trying to estimate a large latent variable model using RStan. I'm getting several divergent transitions after warmup when performing estimation on simulated data and am trying to diagnose the problem using pairs() plots. 

I've looked at scatterplots between each parameter and lp__ using the pairs() function. Most divergent transitions occur below the diagonal, suggesting that the model needs to be reparameterized. But how can I tell which parameters are the source of sampling problems? None of the parameters appear significantly correlated with lp__.

[Bob Carpenter]

The pairs plot by default sets up two chains above the diagonal
and two chains below the diagonal. So all this is saying is that
they're happening on a subset of the chains. You want to see where
the divergences are at to see if you can learn anything.

The ideal pairs plot will look like a bunch of bulls-eyes of
independent unit-scaled posteriors. Anything that smacks of high
correlation or varying scales can be a problem for adaptation or
sampling.

- Bob

> On Feb 20, 2017, at 12:57 PM, David <djtma...@gmail.com> wrote:
>
>
> I'm trying to estimate a large latent variable model using RStan. I'm getting several divergent transitions after warmup when performing estimation on simulated data and am trying to diagnose the problem using pairs() plots.
>
> I've looked at scatterplots between each parameter and lp__ using the pairs() function. Most divergent transitions occur below the diagonal, suggesting that the model needs to be reparameterized. But how can I tell which parameters are the source of sampling problems? None of the parameters appear significantly correlated with lp__.
>

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

[Luis Damiano]

On Monday, 20 February 2017 15:16:42 UTC-3, Bob Carpenter wrote:

The ideal pairs plot will look like a bunch of bulls-eyes of
independent unit-scaled posteriors.  Anything that smacks of high
correlation or varying scales can be a problem for adaptation or
sampling.

- Bob

This small bit of experience should be in the manual. Certainly, it's a technical manual and not a handbook, but that's a very, very intuitive rule to guide the modelling process: try to make all the parameters standarized and uncorrelated.

Maybe too obvios, but it's worth noting this is true for actual parameters to be estimated and not for transformed parameters, right?

[David]

Hi Bob,

Thanks very much for your response. A few questions:

1) What do you mean by "you want to see where the divergences are at"? I'm not sure how to operationalize that phrase.

2) I'm using the product of a gamma with a normal to generate a t. How do I think about comparing the "scales" of these different distributions? Is it better to draw uniforms and transform via inverse cdfs? 

3) Suppose that my prior distributions are independent and I'm using the QR decomposition on covariate matrices. What else could I do to eliminate correlation in the posterior?

Thanks,

David

[Bob Carpenter]

> On Feb 20, 2017, at 2:33 PM, Luis Damiano <damia...@gmail.com> wrote:
>
>
> On Monday, 20 February 2017 15:16:42 UTC-3, Bob Carpenter wrote:
> The ideal pairs plot will look like a bunch of bulls-eyes of
> independent unit-scaled posteriors. Anything that smacks of high
> correlation or varying scales can be a problem for adaptation or
> sampling.
>
> - Bob
>
> This small bit of experience should be in the manual.

Agreed. This will be a key principle in our methodology guides.
We've been learning a lot of this ourselves as we go. Michael's
about to roll out some really great case studies later this week
or next week.

> Certainly, it's a technical manual and not a handbook, but that's a very, very intuitive rule to guide the modelling process: try to make all the parameters standarized and uncorrelated.

The Stan manual's going to be split into two parts:

1. the reference manual with the language spec and the functions doc

2. a revised user's guide that we'll publish with CRC and online
in HTML using bookdown (just like Hadley Wickham's Advanced R book)
(and also as a freely downloadable PDF---we're committed to keeping the
doc open)

> Maybe too obvios, but it's worth noting this is true for actual parameters to be estimated and not for transformed parameters, right?

Correct. Keep in mind that's on the unconstrained scale,
not the constrained scale. So if you have a parameter with
<lower = 0>, it's the log of the parameter. If you have
a covariance matrix, it's the Cholesky factor with log-transformed
diagonal. Unfortunately, you can't access that directly through
RStan because it can't get the transforms out from the language (yet).

- Bob

[Bob Carpenter]

> On Feb 20, 2017, at 2:52 PM, David <djtma...@gmail.com> wrote:
>
> Hi Bob,
>

> Thanks very much for your response. A couple questions:

>
> 1) What do you mean by "you want to see where the divergences are at"? I'm not sure how to operationalize that phrase.

The points at which divergent transitions arise are
highlighted in the plots. Or you can get them directly
in RStan using this:

https://cran.r-project.org/web/packages/rstan/vignettes/stanfit-objects.html

> 2) I'm using the product of a gamma with a normal to generate a t. How do I think about comparing the "scales" of these different distributions? Is it better to draw uniforms and transform via inverse cdfs?

Why do that rather than have the Student-t itself? Multiplicative
paramters are bad news as they tend to lead to banana-like posteriors,
which are terrible --- no good adaptation parameters given that
curvature varies around the parameter.

> 3) Suppose that my prior distributions are independent and I'm using the QR decomposition on covariate matrices. What else could I do to eliminate correlation in the posterior?

That's a question for Ben---you may want to pull it out in
its own message to get the right people's attention.

- Bob

[David]

On Monday, February 20, 2017 at 1:59:59 PM UTC-6, Bob Carpenter wrote:

> On Feb 20, 2017, at 2:52 PM, David <djtma...@gmail.com> wrote:
>
> Hi Bob,
>
> Thanks very much for your response. A couple questions:
>
> 1) What do you mean by "you want to see where the divergences are at"? I'm not sure how to operationalize that phrase.

The points at which divergent transitions arise are
highlighted in the plots.  Or you can get them directly
in RStan using this:

https://cran.r-project.org/web/packages/rstan/vignettes/stanfit-objects.html

So I can "see" the divergences on the pairs() plots (or can extract them directly), but I don't really know what to do with this information. What patterns am I looking for?

> 2) I'm using the product of a gamma with a normal to generate a t. How do I think about comparing the "scales" of these different distributions? Is it better to draw uniforms and transform via inverse cdfs?

Why do that rather than have the Student-t itself?  Multiplicative
paramters are bad news as they tend to lead to banana-like posteriors,
which are terrible --- no good adaptation parameters given that
curvature varies around the parameter. 

 

Section 26.6 in the manual suggests this parameterization in order to avoid the heavy-tailedness of the t-distribution.

 

> 3) Suppose that my prior distributions are independent and I'm using the QR decomposition on covariate matrices. What else could I do to eliminate correlation in the posterior?

That's a question for Ben---you may want to pull it out in
its own message to get the right people's attention.

Ok if I can't figure it out I'll make a separate post.

Thanks!

David

- Bob

[Ben Goodrich]

On Monday, February 20, 2017 at 1:16:42 PM UTC-5, Bob Carpenter wrote:

The pairs plot by default sets up two chains above the diagonal
and two chains below the diagonal.  So all this is saying is that
they're happening on a subset of the chains.  You want to see where
the divergences are at to see if you can learn anything.

This hasn't been true for a couple of years. By default, the draws (across all chains) with below median accept_stat__ are below the diagonal and those with above median accept_stat__ are above the diagonal. If there are divergent transitions above the diagonal, those can usually be fixed with a higher adapt_delta. If they are below the diagonal, you can try increasing adapt_delta but oftentimes you have to adjust the model.

Ben

[Bob Carpenter]

Sorry about that---I hadn't realized the behavior changed.

Also forgot that Ben (I think?) put that reparameterization
in the manual. I haven't tried it, but assume whoever put
it in the manual made sure it worked first!

- Bob

[David]

Hi Ben,

Thanks for your response. I was wondering if there was other information I was suppose to glean from the pairs() plots. In particular, if the divergent transitions are below the diagonal how do I know what part of the model I need to adjust?

Thanks,

David

[Ben Goodrich]

On Monday, February 20, 2017 at 4:44:46 PM UTC-5, David wrote:

Thanks for your response. I was wondering if there was other information I was suppose to glean from the pairs() plots. In particular, if the divergent transitions are below the diagonal how do I know what part of the model I need to adjust?

You want to look below the diagonal for weird shapes with non-constant variance, particularly funnels and the like. This often happens with variance parameters or hyperparameters. Sometimes specifying log = TRUE helps when some of the variables are strictly positive but get very close to zero.

Ben

[andre.p...@googlemail.com]

Hi David,

I often look at n_eff samples in the summary report. I reparameterize where the n_eff is low or change the distribution, from cauchy to student-t with nu = 4.

Sometimes I use the inv_gamma distribution (but I fit the range values with fitdistrplus) R-package and give some extra to the variance.

Sometimes applying the MATT trick does help too. I never figured out however, why there is this problem in my models.

Thus I felt modeling is an art not only a science.

Andre

[Luis Damiano]

> Certainly, it's a technical manual and not a handbook, but that's a very, very intuitive rule to guide the modelling process: try to make all the parameters standarized and uncorrelated.

The Stan manual's going to be split into two parts:

1.  the reference manual with the language spec and the functions doc

2.  a revised user's guide that we'll publish with CRC and online
    in HTML using bookdown (just like Hadley Wickham's Advanced R book)
    (and also as a freely downloadable PDF---we're committed to keeping the
    doc open)

How can I contribute to #2? I'd *really* love to write and plot useful stuff that will help statisticians in their journeys (something applied comprising both computations and statistical data analysis).

[Bob Carpenter]

There are two different things going on here.

For computational efficiency and numerical stability, we'd like
the parameters to be independent in the posterior and unit scale.
Ideally without fat tails. Unit normal works well.

Of course, we can't get that, but often we can reparameterize to
get the same inferences with a better behaved posterior. The centered
parameterization is better when there's a lot of data and the non-centered
when there's less; better here is looking more like a standard normal.

The other issue is the statistical aspect of the model. Often we wind
up with posteriors that are inconsistent with our prior information. Often
we're vague about prior information, but there's enough data to make up
for it. When there's less data, there's much more prior sensitivity.
So we want to look at the posteriors to see if they make sense, and if
not, apply stronger priors so that the end-to-end inferences make sense.

Andrew's "folk theorem" ties these things together---misspecified models
(either in the prior or in the likelihood) are often also hard to fit.

- Bob

[Bob Carpenter]

If I get my way, we'll be developing (2) online in a GitHub
repo. I'm not sure we're going to be able to coordinate
too many authors and we're going to need to come up with
a presentation style that all of us can live with (which will
probably be the hardest part of the whole thing---a classic
bike shed issue, though we'll all deny it).

- Bob

[andre.p...@googlemail.com]

Bob,

ok, does the pairs-plot does it give some valuable information, about the location of the

divergences? Otherwise, we could just say, ok we have divergence, so we increase the adapt 

delta or change the distribution slightly.

When I discovered divergence, I cannot specify a special region of parameter problem

usually.

Andre

[Bob Carpenter]

It plots them in a separate color, I believe, but I'm
clearly not up to date on the pairs plot any more. One of
the plots does this. Then you can see if they're all in
the neck of a funnel or on one side of a banana-shaped
posterior, you can detect some problems.

- Bob

[Jonah Gabry]

On Friday, February 24, 2017 at 9:03:41 PM UTC-5, Bob Carpenter wrote:

It plots them in a separate color, I believe, but I'm
clearly not up to date on the pairs plot any more.  One of
the plots does this.  Then you can see if they're all in
the neck of a funnel or on one side of a banana-shaped
posterior, you can detect some problems.

- Bob 

Yeah divergences are plotted in red and hitting max-treedepth is indicated with yellow. There is also a new ggplot2-based version of the pairs plot that I pretty much have ready to go (https://github.com/stan-dev/bayesplot/pull/69). It will be released as part of the bayesplot package. The ggplot version will also give you more control over the appearance of the points indicating divergences and hitting max-treedepth (sometimes in the base R version they are hard to see).

On another note, the trace plots in the bayesplot package (the mcmc_trace function) allow you to add a rug plot that shows where the divergences occurred. This isn't always helpful, but it gives you a different way to look at the divergences that can complement what the pairs plot tells you. Here's a quick demonstration using rstanarm, but it works the same for models fit using rstan too: 

library("rstanarm")

fit <- stan_glm(mpg ~ ., data = mtcars,
  # next line to keep example fast and also ensure we get some divergences
                prior = hs(), iter = 400, adapt_delta = 0.8)

# extract draws using as.array (instead of as.matrix) to keep
# chains separate for traceplot
posterior <- as.array(fit)

library("bayesplot")

# grab divergences ('nuts_params' is a convenience function in bayesplot)
divs <- nuts_params(fit, pars = "divergent__")
mcmc_trace(posterior, pars = "sigma", divergences = divs)

This results in the plots in the attached png file.  

Jonah

[Jonah Gabry]

Actually to get it to look like plot in the png I attached you'd also need to change the default color scheme by first calling 

bayesplot::color_scheme_set("mix-brightblue-gray")

but that's not essential.