Results difference with different parametrisation: The case of the Anova I random model

[Pierre Lebrun]

Hello all,

I am writing this question to have a discussion on very different results I obtained using two different Stan parametrisation of the same model. I would have thought that my STAN models provide the same results, but unfortunately no. I am also comparing the results obtained with a similar parametrisation with JAGS (JAGS is not dead, right ? :-) ) and with beta-expectation tolerance intervals obtained analytically, using Robert Mee 1984 paper (Mee R.W., Beta-expectation and beta-content tolerance limits for balanced one-way ANOVA random model. Technometrics. 1984;26:251-254.). Notice I am not taking the results of Mee as granted (there are some approximation), but I know by experience that they are a good approximation of the tolerance intervals.

A little bit of background: the target is to obtain the posterior distribution of the parameters and the beta-expectation tolerance interval using MonteCarlo estimation (which, in this case, are the prediction interval obtained by integrating posterior parameters into the model/likelihood). Under flat prior, I am (well,... I was) confident that the posterior modes must be similar to the maximum likelihood estimators, and the beta-expectation tolerance interval of Mee's paper is the same than the MC-computed one.

Next, the data set is one of these that provides difficulties to maximum likelihood estimator, with a negative variance component for the random effect (and thus, lme() or SAS/PROC MIXED would put the estimate to 0).

Here is the results with the first part of the code (attached), really just translating the ANOVA 1 random model into Stan (Note also I obtained  similar results when running the ANOVA models obtained with brms or rstanarm package, for what's worth. I am running Stan 2.9.0-3 on Mac OS X 10.10.5 and R 3.2.2, but also, similar results have been obtained under Windows 7 (Stan 2.9.0 on R 3.2.2).

Observation : posterior mode are well achieved. Variability of b0 (the fixed intercept) is way too big. Prediction interval (beta-exp. tolerance interval) is way too large (right figure, red). Actually the HPD on b0 is already larger than the tolerance interval obtained using Mee's paper (blue).

And here are the results obtained with the parametrisation in terms of precision (as one would do with JAGS, for instance). Note that JAGS chains and posterior densities looks better, but predictive interval is similar using JAGS.

 

Observation: Chains and densities look not good except with stronger prior (which also give better running time). But interval is more plausible (In the sense that I would buy it, if I could not see all the other results). Notice that similar results are obtained with JAGS (model attached at the end of the script, but not the running environnement (what we have is sufficient for the discussion and that's already big I think).

So, my problem is that given similar data, I would except different coding to provide different effective sample size, different qualities of results, but not "different" results... And of course I don't know what is preferable, but I have "frequentist" tools, well-know, used and verified that tell me that none of the solution I obtained is good.

So, the code is fully standalone, so please tell me if you observe similar results, if you see that I am making a stupid coding mistake, or if this should be seen as a bug or something to improve. I admit I tried to run it with older version of Stan but I did not succeed to install them properly so I gave up. Can someone run that under rstan 2.4 or something ?

Thanks in advance !

Pierre

[Bob Carpenter]

You want to simulate data from the model and test that
the model is calibrated in that the expected number of
posterior intervals contain the true values.

And yes, two different parameterizations will give you
the same answer. If they don't, then there's probably an issue
with the transform such as a missing Jacobian adjustment.

- Bob

> On Mar 24, 2016, at 10:55 AM, Pierre Lebrun <lebrun...@gmail.com> wrote:
>
> Hello all,
>
> I am writing this question to have a discussion on very different results I obtained using two different Stan parametrisation of the same model. I would have thought that my STAN models provide the same results, but unfortunately no. I am also comparing the results obtained with a similar parametrisation with JAGS (JAGS is not dead, right ? :-) ) and with beta-expectation tolerance intervals obtained analytically, using Robert Mee 1984 paper (Mee R.W., Beta-expectation and beta-content tolerance limits for balanced one-way ANOVA random model. Technometrics. 1984;26:251-254.). Notice I am not taking the results of Mee as granted (there are some approximation), but I know by experience that they are a good approximation of the tolerance intervals.
>
> A little bit of background: the target is to obtain the posterior distribution of the parameters and the beta-expectation tolerance interval using MonteCarlo estimation (which, in this case, are the prediction interval obtained by integrating posterior parameters into the model/likelihood). Under flat prior, I am (well,... I was) confident that the posterior modes must be similar to the maximum likelihood estimators, and the beta-expectation tolerance interval of Mee's paper is the same than the MC-computed one.
>
> Next, the data set is one of these that provides difficulties to maximum likelihood estimator, with a negative variance component for the random effect (and thus, lme() or SAS/PROC MIXED would put the estimate to 0).
>
> Here is the results with the first part of the code (attached), really just translating the ANOVA 1 random model into Stan (Note also I obtained similar results when running the ANOVA models obtained with brms or rstanarm package, for what's worth. I am running Stan 2.9.0-3 on Mac OS X 10.10.5 and R 3.2.2, but also, similar results have been obtained under Windows 7 (Stan 2.9.0 on R 3.2.2).
>
>
>
>
>
>

> Observation : posterior mode are well achieved. Variability of b0 (the fixed intercept) is way too big. Prediction interval (beta-exp. tolerance interval) is way too large (right figure, red). Actually the HPD on b0 is already larger than the tolerance interval obtained using Mee's paper (blue).
>
> And here are the results obtained with the parametrisation in terms of precision (as one would do with JAGS, for instance). Note that JAGS chains and posterior densities looks better, but predictive interval is similar using JAGS.
>
>
>
>
>

> Observation: Chains and densities look not good except with stronger prior (which also give better running time). But interval is more plausible (In the sense that I would buy it, if I could not see all the other results). Notice that similar results are obtained with JAGS (model attached at the end of the script, but not the running environnement (what we have is sufficient for the discussion and that's already big I think).
>
> So, my problem is that given similar data, I would except different coding to provide different effective sample size, different qualities of results, but not "different" results... And of course I don't know what is preferable, but I have "frequentist" tools, well-know, used and verified that tell me that none of the solution I obtained is good.
>
> So, the code is fully standalone, so please tell me if you observe similar results, if you see that I am making a stupid coding mistake, or if this should be seen as a bug or something to improve. I admit I tried to run it with older version of Stan but I did not succeed to install them properly so I gave up. Can someone run that under rstan 2.4 or something ?
>
> Thanks in advance !
>
> Pierre
>

> --
> 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.
> <for stan group .R>

[Michael Betancourt]

Looking at those trace plots (i.e. where the chains get “stuck” for long periods)
I wouldn’t trust the results from either fit. Something is either wrong with your
model or you need to reparameterize it to remove pathologies.

[Jonah Gabry]

In addition to what Bob and Michael said, a few other comments:

1) Why no prior on b0? 

2) In the model parameterized in terms of the precision try something like cauchy(0, small number) instead of uniform(0, big number), where in your case "small number" isn't tiny but doesn't have to be huge. Maybe 10 or so. I'd prefer to not work with the precision at all, but if you want to then I would avoid the hard constraint of the uniform. 

3) You can speed up your code by vectorizing the likelihood. So in the model block you'd have something like 

model{
  vector[N] tmp;

  // ...

  for (n in 1:N) {
   tmp[n] <-  beta0[Series[n]];
  }
  y ~ normal(b0 + tmp, sig_e);  // vectorized
}

4) You definitely shouldn't need 62500 iterations with 50000 of them warmup, and also I'd say get rid of thin=5. For JAGS these are typical settings but for Stan that's almost always overkill.

Jonah

[Pierre Lebrun]

Hi Bob, Michael and Jonah,

Thanks for your first comments and your speed of reaction ! Very much appreciated ! 

@Bob, I did not look at the coverage for the Bayesian model, but I did for the Mee's interval, they are quite good (but results' quality of Mee decreases when degrees of freedom are low... hence the idea to replace it by a Bayesian model). Testing coverage with simulated data and the Bayesian sampler would take a lot of time, and I can see from the interval that the coverage will not be good.. so, no simulation needed at this point :-) ... For the correction with the Jacobian, I am not sure how to handle that. It seems to me that my two models are identical (the "normal" distribution statements are the same). I just add a transformation for the prior distribution to avoid the sampler to draw arbitrary low values for the SD (that become arbitrary high in term of precision). Anyway, precision values are back-transformed into SD before included in the model. 

@Michael, that's the idea of what I am doing: assessing different types of parametrisation (and that's why I also generate code form brms and look at rstanarm available solutions, for this model). So far I am not lucky because results are too much different.

@Jonah. Thanks for the "speed up" adjustment of the code. Very effective indeed :-) ! Now, for the prior, when I use a brms or rstnarm model, these types of priors are already included, and results were like the "SD" model. I can indeed also improve the "precision" model with the cauchy prior, but I obtain satisfying results only with a cauchy(0,very_big_number)... like 10000000 (so, as the uniform). Otherwise, the posteriors are too biaised and don't match the ML estimates at all. Note, finally, that for the iterations, I know it is overkill, but I increased them (and the thin) because of the problems I have. 

So, a priori my idea was to use no prior at all, and let the data speaks as much as possible for a more formal comparison of Frequentist vs. Bayesian results. I know using no prior is generally not recommended. Really I have a concern if my models are indeed identical. If they are not, I would be delighted to now what can be done to have them identical. 

Anyway, I would generally only use the model parametrised in term of SD. But now I am less confident because in general I do that to obtain the predictive distribution, compute risk to be out of some boundaries, etc. If the predictive is not good then I am in trouble :-)

To recap, I also know that the data set is a very bad one (see note on poor variance REML estimates in first message). But that's the type of data we encounter all the time (that's actually blinded real data). So I don't know what can be done. I have the feeling that previous version of Stan did not suffer from the problem (but that's just a feeling and I cannot check it).

Kind regards and thank you again !

Pierre

 

[Bob Carpenter]

I took a look at the models this time. Is it these two that you
think should give you the same result?

data {
int J;
int N;
real y[N];
int <lower=1> Series[N];
}
parameters {
real b0;
real beta0[J];
real <lower=0> sd_b;
real <lower=0> sd_e;
}
model {
for(i in 1:N)
y[i] ~ normal(b0 + beta0[Series[i]], sd_e);
beta0 ~ normal(0, sd_b);
}

===========================

data {
int J;
int N;
real y[N];
int <lower=1> Series[N];
}
parameters {
real b0;
real beta0[J];
real <lower=0> tau_b;
real <lower=0> tau_e;
}
transformed parameters {
real <lower=0> sig_b;
real <lower=0> sig_e;
sig_b <- sqrt(1/tau_b);
sig_e <- sqrt(1/tau_e);
}
model {
tau_b ~uniform(0,1000000000000);
for(i in 1:N)
y[i] ~ normal(b0 + beta0[Series[i]], sig_e);
beta0 ~ normal(0, sig_b);
}

=======================

The first model provides a uniform improper prior on
scale and the second a mix of improper and fat uniform
on precision. I'm not sure how to think about the change
of variables with improper or narrow priors like this.
uniform(0, a) on scale, though, is not the same as
uniform(0, pow(a, -2)) on precision --- you need the Jacobian
adjustment.

I don't think a big fat uniform is going to do much in practice.
In general, Stan models should have finite log density for any
parameter values satisfying the constraints, which your second model
doesn't do (probably won't matter here if the posterior's proper).

You can just use improper uniform priors in Stan, as you
can see for tau_e. I don't know what you thought the fat uniform
on tau_b would do.

Posterior means only approach the MLE estimates in the limit.
If you want MLE estimates, you can just do MLE. But these hierarchical
models often have no MLE, with growing likelihood as sig_b -> 0 and
beta0 -> 0.

There hasn't been any change in Stan's behavior as far as sampling
goes for a couple years, and the change was added robustness.

- Bob

[Pierre Lebrun]

Thanks Bob!

As you mentioned, we have this behaviour of the likelihood as sig_b -> 0, and MLE actually also fails by "resetting" the value at sig_b=0 after obtaining negative values. The reason why I am not using MLE is more general: I want to obtain predictive intervals. For simple model like the random ANOVA I model, I can have nice approximations (the Mee's paper I mentioned), but when designs are sparse and unbalanced (Pharma industry problem.. occurs often) these can become tricky. Let apart a more complex mixed-effect models. There, the Bayesian way is the only one.

So, I wanted to obtain the same results, but of course, without the prior statement of sig_b in the second model.

For the prior, initally, I did not want to use any prior on sig_b (I firstly used the default improper prior)... but without it, Stan failed to converge for the second model  (this is also the reason why I use absurdly high burnin and iterations). Thus, I tried constraining the uniform with trial and error (largest value so the model converge).. I guess the solution is to use the first model with reasonable priors (that I have, in this case, but I would have liked to be more general and stick to improper priors)... also, due to that, I should not look so formally to the coverage I can get (because coverage is dependant of my prior in this case).

So just to rephrase my concern: I want to use the first model (the second model was to mimic a JAGS model). It actually runs without a convergency issue, all model checks are on the green (except I can see the uncertainty of the posteriors being too large). But the predictive interval is not of beta-expectation as I would have imagined, and to a big margin ! (notice that, with better data sets, it is of beta-expectation !). Here, I can discover the problem because I can cross-check using Mee's paper and MLE... but with other more complex cases, I would provide too large intervals without easily seeing they do not predict well the future uncertainty.

So, thanks anyway, conclusion is that we'd better be very careful with such poor data :-)

Happy Easter to you and your team !

Pierre

[Bob Carpenter]

Our advice on priors is in the manual chapter on regression
and here:

https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations

Poor data should just mean wide posterior
intervals --- the model should still be well-calibrated for
data generated (roughly) according to the model's assumptions.

- Bob