Log-normal parameter: can I use the mean in stan object?

Modeling
1

Hi,

The question:

In my hierarchical model the parameters of interest are log-normally distributed, so do I still use the ‘mean’ estimates provided in the stan object after fitting as my estimate for the parameter - that is, the ‘mean’ estimates you see when running ‘print(object)’? Or am I best to calculate the mode of the (log-normal) distribution of all my draws for each parameter?

In case it is useful, here is the code snipped defining the parameters as log-normal; the full code is included at the end of the post.

noise ~ lognormal(meannoise, sdnoise);

  alpha ~ lognormal(meanalpha, sdalpha);

  M ~ lognormal(meanM, sdM);

The parameter of interest is ‘M’ and it is estimated separately for each individual completing my experiment.

BTW, I previously posted my model and received helpful assistance regarding ‘divergence’ before: RStan divergence issue - assistance with experimental research - I am happy to say ‘divergences’ are not a problem any longer: thank you!

Thanks for having a look and let me know if you have any further questions that I can clarify.

Cheers,

Alex

****

Here’s the full Stan model code:

data {

  int N; // number of trials total (across participants), integer

  int nsubj;  // number of subjects,

  int choices[N]; // the binary choice vector

  real <lower=0> x2a[N];

  real <lower=0> x1b[N];

  real <lower=0> p2a[N];

  real <lower=0> p1b[N];

  int sid[N];

}



parameters {

  // Group-level:

  real<lower=-2.30, upper=1.61> meannoise; // assuming mean is 0.1 to 5 -  

  real<lower=0, upper=1.13> sdnoise; // according sds  - calculated as sd = (b - a) / sqrt(12)

  real<lower=-2.3, upper=2.3> meanalpha; // assuming mean is 0.1 to 10

  real<lower=0, upper=1.33> sdalpha;

  real<lower=-2.3, upper=4.5> meanM; // assuming mean is 0.1 to 90

  real<lower=0, upper=1.96> sdM;

  

  // Individual-level:

  real<lower=0.1> noise[nsubj]; // Noise, constrained it to be above [0.1, INF]

  real<lower=0.1> alpha[nsubj]; // alpha, constrained it to be above [0.1, INF]

  real<lower=0.1> M[nsubj];  // M, constrained it to be [0.1, INF)

}



model {

  real ua; // utility of the option a

  real ub; // utility of the option b

  // Group-level parameters:

  meannoise ~ uniform(-2.30,1.61); // assuming mean is 0.1 to 5

  sdnoise ~ uniform(0,1.13); // according sds

  meanalpha ~ uniform(-2.3,2.3); // assuming mean is 0.1 to 10

  sdalpha ~ uniform(0,1.33);

  meanM ~ uniform(-2.3, 4.5); // assuming mean is 0.1 to 90

  sdM ~ uniform(0,1.96);

  

  // Individual-level parameters:

  noise ~ lognormal(meannoise, sdnoise);

  alpha ~ lognormal(meanalpha, sdalpha);

  M ~ lognormal(meanM, sdM);

    

  for (i in 1:N) {

      int t = sid[i];

      ua=0.0;

      ua += p2a[i]*((x2a[i]^alpha[t])/((x2a[i]^alpha[t])+(M[t]^alpha[t])));

      

      ub=0.0;

      ub += p1b[i]*((x1b[i]^alpha[t])/((x1b[i]^alpha[t])+(M[t]^alpha[t])));



      choices[i] ~ bernoulli_logit(((ua-ub)*noise[t]));

      }

}
2

Mean, median, mode–the world is your oyster. Pick the one you want. I’ve seen all three recommended in various Bayesian textbooks. If you’re worried about it, report two or even all three. You could also display the entire posterior distribution of your parameter[s] of interest.

2 Likes
3

Are you noticing that in practice the posterior looks lognormal, and want a good way to summarize the posterior?

Or are you noticing that the hierarchical prior on the parameter M is lognormal, and want a good way to summarize the prior distribution?

Or are you simply noticing that the prior on M is lognormal, and generally worried about whether this causes problems in reporting?

1 Like
4

Thank you Solomon. I have done some more reading and I’m converging on the same opinion that it kind of depends and all are valid, and perhaps best to report all three.

1 Like
5

It’s actually not the prior M that I’m that interested in, or what I call the ‘group-level’ parameter M; instead it’s the parameter M for each individual (‘individual-level’). This latter parameter is log-normally distributed by-design, and I was wondering if I should use the ‘mean’ provided in the output by RStan or calculate the moment (e.g., calculate the modes or medians of all the draws) myself when it’s log-normal.

6

The log-normal distribution by design is over the population of parameters M. The posteriors for each individual element of M are not necessarily lognormal. You can summarize them in whatever way makes sense, as @Solomon says. Often such a summary would include quantiles or a credible interval. But there’s nothing particularly special about the fact that the prior is lognormal (really a mixture of lognormals) that changes the strategy for reporting these margins of the posterior.

1 Like
7

Thank you jsocolar for explaining: that’s helpful point about the individual versus population of parameters. I will make sure to report the margins of the posterior in my summary as this will make it clearer.

8

A lot depends on how you or others will be using the posterior and its summaries.

In my field (cosmology) we have a set of interesting parameters (power spectra) which are positive definite. This means that their posteriors are always skew positive, and indeed look a bit like a lognormal. In many cases these distributions are, nonetheless, very well approximated as multivariate normal distributions (with a cutoff at 0), but centred at the mode of the posterior.

(Unfortunately, the mode is very difficult to find in high-dimensional problems, as the typical set avoids it due to the measure factor.)

9

Thanks, your comment prompted me to look into the log distribution a bit more and seeing its quite common in nature.