Random-variable objects; functions to grab posterior means, medians, sds, and mad_sd's from stan objects; also more general versions of these ideas

[Andrew Gelman]

Hi, I have a few questions/requests which are relevant for my applied work, and also for the new book with Jennifer.

1. Random-variable objects.

After I've fit a model, I'll often want to do things with the parameters. For example, MRP. I do it in R. But it's awkward to work with posterior simulations in R for all the usual reasons: I'm continually having to wrestle with array algebra, having to decide what to do with indexes. For example, if I want to draw from normal(X*beta, sigma), I need to either figure out some slick "apply" thing or else loop on simulation draws and use indexing like beta[s,] and sigma[s], which also required me first assigning a blank object to whatever I'm creating. Yes, it can be done, but it turns simple statistical manipulations, which would be short one-liners with random variable objects, into mini-programs. This is, I think, contrary to the spirit both of R and of probabilistic programming.

In Stan, we can work with random variables in the generated quantities block. But right now you can't run a generated quantities block on its own. I think this should change: we should allow the user to pass a stanfit object (or whatever this is called outside of R, umm, I guess it's just "the output from a fitted Stan model") and run a generated quantities block all on its own.

I suppose the alternative would be to do this using Stan functions but then we'd need a workflow for that. Right now I find it more difficult to write a function than to write a generated quantities block, but maybe that would change with practice.

The real point is that I'd like to do math operations with random variables, and Stan seems like the right language to do this in. Ten years ago, Jouni programmed random-variable objects in R, but as of now that does not seem like the right way to go, because (a) R is clunky, and (b) now we already have Stan which does this.

From a Stan perspective, this is also good stuff because it's making use of Stan's capability of being a programming language (rather than the usual thing, which is that people think of Stan as a graphical modeling language, even though it isn't). Anyway, this would be a big help both for applied work and for teaching / textbook writing.

2. Stan output.

Currently, in addition to Rhat and n_eff, we report posterior mean, sd, and a bunch of quantiles. I'm wondering if we should follow the model of rstanarm and replace mean and sd by median and mad_sd (that's the median absolute deviation from the median, multiplied by a scaling factor so that it is equivalent to the sd for the case of the normal distribution). These have the advantage that they always exist (so will work, for example, with the Cauchy distribution), also I think they're closer to what users should really want in most applications. Similarly, we'd replace se_mean by mad_se which is just mad_sd/sqrt(n_eff). Then for the quantiles, now we report 2.5%, 25%, 50%, 75%, 97.5%. OK, we don't need 50% anymore. I'd recommend removing 2.5% and 97.5% too (or maybe not, I'm not quite sure; on one hand, focusing on these intervals is bad practice, and these tail values have lots of simulation noise; on the other, it can be good for people to see some wide uncertainty bounds).

So here's the new output:

median mad_sd 25% 75% n_eff mad_se Rhat

You can see that I changed the order a bit. I find se_mean (or, equivalently, mad_se) a bit difficult to explain. I end up just saying it's sd/sqrt(n_eff) (or mad_sd/sqrt(n_eff), but, if so, it makes sense for it to appear _after_ n_eff in the table. Finally, I kept Rhat at the end because it's convenient to scan the output and look at R_hat.

3. Pulling out parameter estimates:

Suppose I want to just grab point estimates and se's from a stanfit object in R. I can do something like this:
beta_hat <- summary(fit, pars="beta")$summary$mean
(Actually I'd prefer median but that's no big deal right now.)
But the above is quite a mouthful.
Also, if beta is a matrix or array of 2 or more dimensions, I'd like to be able to automatically grab the estimates in matrix or array form, which the above won't do for me.

Or I could do it directly from the Stan output, but that involves looping or figuring out apply() and also doesn't easily generalize beyond scalars and 1-way arrays.

I _think_ what I'm looking for is a function such as posterior_mean(fit, pars="beta")) which returns an object the same size and shape as beta, or maybe something like posterior_mean(fit, pars=c("alpha,"beta","gamma")) which returns a list.

Another way would be something like this: posterior_summary(fit, pars, f), which would apply the function "f" to the parameter or parameters in "pars". For example, posterior_summary(fit, "beta", median) or posterior_summary(fit, c("alpha","beta"), mad).

Hmmm, now I'm thinking that this idea could get us in the direction of working with random variables! For example, what about this sort of function:

posterior_computation(fit, f, data_args, rv_args). Here, "f" is a user-supplied function and it is applied to the stanfit object "fit", along with optional arguments "data_args" and "rv_args". The "data_args" are data supplied to the "foo" function; the "rv_args" are also supplied to "foo" but as random variable objects, meaning that they are grabbed one simulation draw at a time). For example:
foo <- function(X, beta, sigma) {y_hat <- X%*%beta; N <- length(y_hat); return(rnorm(N, y_hat, sigma))};
sims <- posterior_computation(fit, foo, data_args=X_new)
If "foo" returns an array of dimensions "dim_foo", then sims will be an array of dimensions (n_sims, dim_foo), where n_sims is the number of simulations in the stanfit object "fit".

The above will not work in complete generality because I was assuming that "foo" would return an array of fixed dimension (so that they can be put together in one big array that's one dimension higher). But in general "foo" could return anything. So I guess that in general, posterior_computation() would have to return a list of length n_sims. Such a list would be more awkward to work with so maybe there could be an option to collapse it back into an array if all the elements in the list are scalars or else arrays of identical size and shape.

Also we'd want these all to work for stanreg objects but that should be easy enough, I assume, as stanreg objects are themselves the output of stan fits, with some extra information.

4. Summary

In rstan, I'd like to see these new functions (or maybe they exist already?):

posterior_summary(fit, pars, function)
posterior_computation(fit, function, data_args, rv_args)
and some sort of collapsing function for random variable lists.

It would also be good to be able to do more of this in Stan, using the generated quantities block (or maybe using functions but then I need a clean workflow for this in R).

Any thoughts?

[Dustin Tran]

Not to self-advertise too much, but I recommend looking at Edward (http://edwardlib.org) regarding (1); see http://edwardlib.org/api/model in particular. It could be a helpful reference if similar design specs appear for Stan. 

Edward's modeling language is built around random variables as core building blocks, and performing computation on top of them (e.g., sample from them, apply math operations, stack them). A model is never represented explicitly: rather, you perform inference over any collection of random variables, meaning you can interplay model/inference/model checking as you'd like. This lets you loop over all steps of BDA in the same language.

Dustin

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

[Ben Goodrich]

On Monday, February 20, 2017 at 7:04:34 AM UTC-5, Andrew Gelman wrote:

Any thoughts?

That is a lot to think about. Some of it was already planned to some extent. Let's talk about it more soon.

Ben

[Mike Lawrence]

I'm not sure this solves all of Andrew's wish-list, but I recently created a few helper functions to ease the workflow for regression-like models at least. 

Attached is the code, which first has a special contrast matrix function that wraps model.matrix but additionally stores the data as an attribute of the resulting matrix, that you would use as:

  X = get_contrast_matrix(

  data = my_data

  , formula = ~ age * sex

  )

Then a function to use this matrix plus the stan_fit to work out posterior predictions on the conditions of the experiment (assumes the model was fit with the matrix form of regression as on pg 121 of the manual):

  conditions = get_condition_post(

    post = post_from_stan

    , par = 'beta'

    , X = X

  )

Which you can then work with in a standard dplyr chain, for example:

  conditions %>%

    group_by(

        age

        , sex

    ) %>%

    summarize(

        med = median(value)

        , lo95 = quantile(value,.025)

        , lo50 = quantile(value,.25)

        , hi50 = quantile(value,.75)

        , lo95 = quantile(value,.975)

    ) 

--
Mike Lawrence
Graduate Student
Department of Psychology & Neuroscience
Dalhousie University

~ Certainty is (possibly) folly ~

[Michael Betancourt]

It’s fine if you want to add MAD, although I don’t find it as useful, but 

we absolutely cannot remove mean and mcmc_se as these are the

only theoretically-meaningful outputs of the Markov chains.  There is,

for example, no CLT for mad_d/sqrt(n_eff).  You could strain and

find something sort of similar but it would have some exceptional

behavior making things awkward.

2.5% and 97.5%, however, should definitely go.

On Feb 20, 2017, at 7:04 AM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

[Bob Carpenter]

> On Feb 20, 2017, at 7:04 AM, Andrew Gelman <Gel...@stat.columbia.edu> wrote:
>
> Hi, I have a few questions/requests which are relevant for my applied work, and also for the new book with Jennifer.
>
> 1. Random-variable objects.
>
> After I've fit a model, I'll often want to do things with the parameters. For example, MRP. I do it in R. But it's awkward to work with posterior simulations in R for all the usual reasons: I'm continually having to wrestle with array algebra, having to decide what to do with indexes. For example, if I want to draw from normal(X*beta, sigma), I need to either figure out some slick "apply" thing or else loop on simulation draws and use indexing like beta[s,] and sigma[s], which also required me first assigning a blank object to whatever I'm creating. Yes, it can be done, but it turns simple statistical manipulations, which would be short one-liners with random variable objects, into mini-programs. This is, I think, contrary to the spirit both of R and of probabilistic programming.

To code this up and R, you'd need to reimplement all the functions you want
to work with random variable data types. This is what that package with Jouni.

> In Stan, we can work with random variables in the generated quantities block. But right now you can't run a generated quantities block on its own. I think this should change: we should allow the user to pass a stanfit object (or whatever this is called outside of R, umm, I guess it's just "the output from a fitted Stan model") and run a generated quantities block all on its own.

Yes. Everyone agrees this is a good feature. We'll add it
pretty soon, but not in 2.15---that's going to be just getting
the refactor into production.

Just need to get the round trip of converting outputs into
initializations.

...

> 2. Stan output.
>
> Currently, in addition to Rhat and n_eff, we report posterior mean, sd, and a bunch of quantiles. I'm wondering if we should follow the model of rstanarm and replace mean and sd by median and mad_sd (that's the median absolute deviation from the median, multiplied by a scaling factor so that it is equivalent to the sd for the case of the normal distribution).

That sounds too complicated to me. I always think of these things
in terms of how we want to document it.

> These have the advantage that they always exist (so will work, for example, with the Cauchy distribution), also I think they're closer to what users should really want in most applications. Similarly, we'd replace se_mean by mad_se which is just mad_sd/sqrt(n_eff). Then for the quantiles, now we report 2.5%, 25%, 50%, 75%, 97.5%. OK, we don't need 50% anymore. I'd recommend removing 2.5% and 97.5% too (or maybe not, I'm not quite sure; on one hand, focusing on these intervals is bad practice, and these tail values have lots of simulation noise; on the other, it can be good for people to see some wide uncertainty bounds).
>
> So here's the new output:
>
> median mad_sd 25% 75% n_eff mad_se Rhat

Are you going to stop using means altogether? I think there'd
be a user revolt if we don't have mean and se_mean (though I'll
admit users get confused between se_mean and posterior sd).

> You can see that I changed the order a bit. I find se_mean (or, equivalently, mad_se) a bit difficult to explain.


> I end up just saying it's sd/sqrt(n_eff)

That seems pretty easy to me.

> (or mad_sd/sqrt(n_eff), but, if so, it makes sense for it to appear _after_ n_eff in the table. Finally, I kept Rhat at the end because it's convenient to scan the output and look at R_hat.

I don't mind the order.

> 3. Pulling out parameter estimates:

I'll leave this up to the RStan devs. I don't have problems with
this the way things are, but can see how it could be streamlined
for some common operations.

- Bob

[Andrew Gelman]

The reason for using the mad_sd is that it gives roughly the "right" behavior for long-tailed distributions such as the t.  It's not a big deal but in writing the new regression book I've liked mad_sd as it's more interpretable for regressions when degrees of freedom are low.