Individual study-level estimates off-centre from population-level estimates in hierarchical model
57 views
Skip to first unread message
Ben Lambert
Aug 24, 2016, 9:25:35 PM8/24/16
to Stan users mailing list
Hello,
I am running a meta-analysis model of mark-release-recapture experiments. I am encountering an issue where the estimates of my parameters from the individual experiments are often skewed away from the centre of my population level estimates.
As an example here is the output from my model applied to data I simulated for 3 experiments,
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
kappa[1] 1.16 0.09 0.96 0.32 0.49 0.81 1.47 3.70 126 1.05
kappa[2] 2.73 0.11 1.48 0.84 1.74 2.41 3.44 6.33 183 1.01
kappa[3] 1.79 0.06 0.95 0.61 1.08 1.60 2.23 4.06 219 1.02
a 0.14 0.00 0.04 0.06 0.11 0.14 0.16 0.21 165 1.02
p 0.79 0.02 0.39 0.26 0.52 0.73 0.99 1.82 276 1.02
r 20.66 0.19 3.33 14.66 18.41 20.33 22.44 28.11 307 1.00
cDecay[1] 0.15 0.01 0.06 0.02 0.11 0.15 0.18 0.28 119 1.05
cDecay[2] 0.16 0.00 0.05 0.06 0.12 0.15 0.19 0.28 259 1.02
cDecay[3] 0.17 0.00 0.06 0.07 0.13 0.16 0.20 0.29 282 1.00
psi[1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 142 1.02
psi[2] 0.03 0.00 0.01 0.01 0.02 0.03 0.03 0.05 274 1.01
psi[3] 0.08 0.00 0.03 0.04 0.06 0.07 0.10 0.14 297 1.00
sigma 0.05 0.00 0.02 0.01 0.03 0.05 0.06 0.09 182 1.02
overallKappa 1.56 0.09 1.74 0.03 0.37 0.94 2.09 6.23 350 1.01
overallcDecay 0.14 0.00 0.06 0.03 0.10 0.14 0.18 0.27 238 1.00
lifeExpectancy[1] 21.05 12.17 118.60 3.62 5.42 6.85 9.38 47.58 95 1.04
lifeExpectancy[2] 7.36 0.26 3.85 3.57 5.28 6.63 8.30 15.81 214 1.02
lifeExpectancy[3] 7.01 0.25 3.78 3.47 4.92 6.22 7.88 14.27 234 1.00
overallPsi 0.05 0.00 0.05 0.00 0.01 0.03 0.07 0.16 400 0.99
overallLifeExpectancy 334.86 320.04 6400.89 3.73 5.63 7.02 9.80 33.28 400 1.00
lp__ 3119.85 0.30 3.02 3112.94 3118.16 3120.19 3122.05 3124.71 99 1.05
and in particular we see that the overall estimates of the parameter 'cDecay' are skewed slightly away from the average effect -- given by 'overallcDecay'.
Whilst the above isn't a particularly strong effect, in my actual data I sometimes see a really big difference. The magnitude of this difference does decrease a bit as I increase the # of simulated experiments in the meta-analysis, although doesn't disappear entirely; in one of my experimental 'batches' I have 20 experiments and the population-level estimates still differ from the bulk of the study-level estimates!
Does anyone have an idea as to why I am seeing the difference between the study-level estimates and the population-level ones?
My model is below. I assume that there is a constant death rate of the organisms, and hence use an exponential survival function. I would have liked to move to a non-centered model, but since most of my parameters are non-negative I couldn't think of a way to do this that worked well (I tried doing the Matt-Trick on lognormal -- rather than normal -- priors, but this didn't work as well as the below.)
data {
int<lower=0> N; ## Total number of observations
int<lower=0> K; ## Total number of groups
int Y[N]; ## Vector of observations for the release
int S[K]; ## Array of ints for the group sizes
vector[N] t; ## Vector of time observations (days)
int Pos[K]; ## Array of the starting position of each dataset
vector[K] RelNumber; ## Vector of the numbers released across each study
}
parameters {
real<lower=0> kappa[K];
real<lower=0> a;
real<lower=0> p;
real<lower=0> r;
real<lower=0> cDecay[K]; ## mortality hazard
real<lower=0,upper=1> psi[K]; ## prob of recapture on a given day
real<lower=0> sigma;
}
model {
for (i in 1:K) {
real lambda[S[i]];
real numberCaught;
numberCaught = 0;
for (j in 1:S[i]){
lambda[j]= psi[i]*(RelNumber[i]-numberCaught)*exp(-cDecay[i]*(t[Pos[i] + j - 1]));
Y[Pos[i] + j - 1] ~ neg_binomial_2(lambda[j]+1e-5,kappa[i]);
numberCaught = numberCaught + Y[Pos[i] + j - 1];
}
}
for (i in 1:K) {
cDecay[i] ~ normal(a,sigma);
kappa[i] ~ exponential(p);
psi[i] ~ exponential(r);
}
a ~ normal(0.09,0.045);
sigma ~ exponential(20);
p ~ lognormal(1,1);
r ~ gamma(40,2);
}
generated quantities {
real overallKappa;
real overallcDecay;
real lifeExpectancy[K];
real overallPsi;
real overallLifeExpectancy;
for(i in 1:K){
lifeExpectancy[i] = 1/cDecay[i];
}
overallcDecay = normal_rng(a,sigma);
while(overallcDecay<0) overallcDecay = normal_rng(a,sigma);
overallKappa = exponential_rng(p);
overallPsi = exponential_rng(r);
overallLifeExpectancy = 1/overallcDecay;
}
Bob Carpenter
Aug 24, 2016, 9:49:17 PM8/24/16
to stan-...@googlegroups.com
When you have truncated parameters and hierarchical models,
you need to include the truncation because the normalizing term
is no longer constant.
parameters:
> real<lower=0> cDecay[K]; ## mortality hazard
> real<lower=0> a;
model:
> cDecay[i] ~ normal(a,sigma);
So that should be
cDecay[i] ~ normal(a, sigma) T[0, ];
Specifically, the log(1 - normal_cdf(0 | a, sigma))
from teh truncation won't be constant---it depends on
parameters a and sigma.
You're also going to have an issue with everything being positive
and having so much posterior uncertainty. You see some hint of
that skew from the difference between the posterior means and medians
(I mean the usual sense of "skew" here, not what you were reporting
about the hierarchical parameters vs. lower-level parameters). In
Bayes you need to integrate, not look at the point estimates.
- 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.
you need to include the truncation because the normalizing term
is no longer constant.
parameters:
> real<lower=0> cDecay[K]; ## mortality hazard
model:
> cDecay[i] ~ normal(a,sigma);
So that should be
cDecay[i] ~ normal(a, sigma) T[0, ];
Specifically, the log(1 - normal_cdf(0 | a, sigma))
from teh truncation won't be constant---it depends on
parameters a and sigma.
You're also going to have an issue with everything being positive
and having so much posterior uncertainty. You see some hint of
that skew from the difference between the posterior means and medians
(I mean the usual sense of "skew" here, not what you were reporting
about the hierarchical parameters vs. lower-level parameters). In
Bayes you need to integrate, not look at the point estimates.
- 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.
Ben Lambert
Aug 25, 2016, 7:52:10 AM8/25/16
to Stan users mailing list
Thanks Bob - have added the truncation now.
I'm still having the same issue, though - the individual study estimates are often nowhere near the population-level ones.
Any ideas on reparameterisations that might avoid this issue?
Best,
Ben
Bob Carpenter
Aug 25, 2016, 8:31:18 AM8/25/16
to stan-...@googlegroups.com
Do you mean the point estimates? Mean or median (50% quantile)?
I was trying to say that this may be the "right" answer given
the skew of positive-only distributions with mass near zero---the
means will inevitably be larger than the modes (if the mode exists).
And as you can see, even the means and medians differ considerably.
The posterior mean is the point estimate that minimizes expected
square error, not the one that is "most likely" in the MLE sense.
I'd try simulating data to see what you get if that's not what
you're doing already.
- Bob
I was trying to say that this may be the "right" answer given
the skew of positive-only distributions with mass near zero---the
means will inevitably be larger than the modes (if the mode exists).
And as you can see, even the means and medians differ considerably.
The posterior mean is the point estimate that minimizes expected
square error, not the one that is "most likely" in the MLE sense.
I'd try simulating data to see what you get if that's not what
you're doing already.
- Bob
Ben Lambert
Aug 25, 2016, 8:57:02 AM8/25/16
to Stan users mailing list
Hi Bob,
Both the mean and median are not aligned -- this is with fake data. With the real data I see a similar thing, but sometimes the differences are quite a lot bigger than what I get in my simulated data... Any ideas?
Best,
Ben
On Thu, Aug 25, 2016 at 1:30 PM, Bob Carpenter <ca...@alias-i.com> wrote:
Do you mean the point estimates? Mean or median (50% quantile)?
I was trying to say that this may be the "right" answer given
the skew of positive-only distributions with mass near zero---the
means will inevitably be larger than the modes (if the mode exists).
And as you can see, even the means and medians differ considerably.
The posterior mean is the point estimate that minimizes expected
square error, not the one that is "most likely" in the MLE sense.
I'd try simulating data to see what you get if that's not what
you're doing already.
- Bob
> On Aug 25, 2016, at 1:52 PM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:
>
> Thanks Bob - have added the truncation now.
>
> I'm still having the same issue, though - the individual study estimates are often nowhere near the population-level ones.
>
> Any ideas on reparameterisations that might avoid this issue?
>
> Best,
>
> Ben
>
> --
> 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+unsubscribe@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/itcTZmFAfFE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
Bob Carpenter
Aug 25, 2016, 9:17:18 AM8/25/16
to stan-...@googlegroups.com
Try playing around with this:
data {
int<lower=0> N;
real<lower=0> y[N];
}
parameters {
real<lower=0> beta; // mean
}
transformed parameters {
real<lower=0> lambda; // parameter
lambda = 1 / beta;
}
model {
y ~ exponential(lambda);
}
with some simulated data (say N = 5; y = runif(N, 0, 1)).
You'll see that the posterior mean won't match the mean
of the data y with MCMC (it's "skewed" in your sense, but the
real skew is coming from the underlying exponential and what
you're calling skew is just a property of the posterior mean
not following the posterior mode). It will match with an MLE.
- Bob
> On Aug 25, 2016, at 2:56 PM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:
>
> Hi Bob,
>
> Both the mean and median are not aligned -- this is with fake data. With the real data I see a similar thing, but sometimes the differences are quite a lot bigger than what I get in my simulated data... Any ideas?
>
> Best,
>
> Ben
>
> On Thu, Aug 25, 2016 at 1:30 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> Do you mean the point estimates? Mean or median (50% quantile)?
> I was trying to say that this may be the "right" answer given
> the skew of positive-only distributions with mass near zero---the
> means will inevitably be larger than the modes (if the mode exists).
> And as you can see, even the means and medians differ considerably.
> The posterior mean is the point estimate that minimizes expected
> square error, not the one that is "most likely" in the MLE sense.
>
> I'd try simulating data to see what you get if that's not what
> you're doing already.
>
> - Bob
>
> > On Aug 25, 2016, at 1:52 PM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:
> >
> > Thanks Bob - have added the truncation now.
> >
> > I'm still having the same issue, though - the individual study estimates are often nowhere near the population-level ones.
> >
> > Any ideas on reparameterisations that might avoid this issue?
> >
> > Best,
> >
> > Ben
> >
> > --
> > 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.
data {
int<lower=0> N;
real<lower=0> y[N];
}
parameters {
real<lower=0> beta; // mean
}
transformed parameters {
real<lower=0> lambda; // parameter
lambda = 1 / beta;
}
model {
y ~ exponential(lambda);
}
with some simulated data (say N = 5; y = runif(N, 0, 1)).
You'll see that the posterior mean won't match the mean
of the data y with MCMC (it's "skewed" in your sense, but the
real skew is coming from the underlying exponential and what
you're calling skew is just a property of the posterior mean
not following the posterior mode). It will match with an MLE.
- Bob
> On Aug 25, 2016, at 2:56 PM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:
>
> Hi Bob,
>
> Both the mean and median are not aligned -- this is with fake data. With the real data I see a similar thing, but sometimes the differences are quite a lot bigger than what I get in my simulated data... Any ideas?
>
> Best,
>
> Ben
>
> On Thu, Aug 25, 2016 at 1:30 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> Do you mean the point estimates? Mean or median (50% quantile)?
> I was trying to say that this may be the "right" answer given
> the skew of positive-only distributions with mass near zero---the
> means will inevitably be larger than the modes (if the mode exists).
> And as you can see, even the means and medians differ considerably.
> The posterior mean is the point estimate that minimizes expected
> square error, not the one that is "most likely" in the MLE sense.
>
> I'd try simulating data to see what you get if that's not what
> you're doing already.
>
> - Bob
>
> > On Aug 25, 2016, at 1:52 PM, 'Ben Lambert' via Stan users mailing list <stan-...@googlegroups.com> wrote:
> >
> > Thanks Bob - have added the truncation now.
> >
> > I'm still having the same issue, though - the individual study estimates are often nowhere near the population-level ones.
> >
> > Any ideas on reparameterisations that might avoid this issue?
> >
> > Best,
> >
> > Ben
> >
> > --
> > You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/itcTZmFAfFE/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/itcTZmFAfFE/unsubscribe.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> 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.
Ben Lambert
Aug 25, 2016, 11:45:14 AM8/25/16
to Stan users mailing list
Hi Bob,
Thanks for your example here. Yes, I can see in my data that some of the population-level posteriors are skewed -- typically towards the prior. This isn't true for all the cases where there is a difference however. See the attached for an example of relatively non-skewed posteriors -- the top is the population-level distribution, and the bottom are two of the individual study estimates. All the individual estimates are comparable to those that I show here, specifically none have their centres near the population-level estimates.
This is an issue for me, because at the moment I can't really justify using the population-level estimates as a summary of the individual experiments. My experimental collaborators are -- quite rightly -- concerned that the population-level estimates aren't anything but a statistical mirage.
The idea of my collaborators is just to use a non-hierarchical model, and just graph the individual estimates. I want to avoid this if possible, because I really see the value in this work as synthesising the data from previous works. However at the moment I think it might be the least bad option.
Any thoughts?
Ben
Bob Carpenter
Aug 25, 2016, 1:41:17 PM8/25/16
to stan-...@googlegroups.com
Yes, don't concentrate on point estimates.
- Bob
> On Aug 25, 2016, at 5:45 PM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:
>
> Hi Bob,
>
> Thanks for your example here. Yes, I can see in my data that some of the population-level posteriors are skewed -- typically towards the prior.
I still don't understand what you mean by "skewed" here. The
term "skew" has a technical meaning in statistics, and I have
a feeling you may be confusing it with the technical term "bias".
> This isn't true for all the cases where there is a difference however. See the attached for an example of relatively non-skewed posteriors -- the top is the population-level distribution, and the bottom are two of the individual study estimates. All the individual estimates are comparable to those that I show here, specifically none have their centres near the population-level estimates.
By "centre" you mean median?
> This is an issue for me, because at the moment I can't really justify using the population-level estimates as a summary of the individual experiments. My experimental collaborators are -- quite rightly -- concerned that the population-level estimates aren't anything but a statistical mirage.
Treating point estimates with certainty is the mirage---you
want to look at the posterior for the whole model. Assuming
that population level point estimates (you haven't said which
point estimates you're using---MLE, posterior median, posterior
mean) are representative of the population is also not right.
See, for example:
https://en.wikipedia.org/wiki/Ecological_fallacy
which in turn cites:
https://en.wikipedia.org/wiki/Simpson%27s_paradox
> The idea of my collaborators is just to use a non-hierarchical model, and just graph the individual estimates. I want to avoid this if possible, because I really see the value in this work as synthesising the data from previous works. However at the moment I think it might be the least bad option.
>
> Any thoughts?
From a practical point of view, if you need to publish with collaborators,
then do the expedient thing for them and the editors for your career and
make it a longer-term goal to get the stats "right". I'm not just being
sarcastic here---it's what I have to do, too. See, e.g.,
https://lingpipe-blog.com/2014/10/29/beckys-and-my-annotation-paper-in-tacl/
- Bob
- Bob
> On Aug 25, 2016, at 5:45 PM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:
>
> Hi Bob,
>
> Thanks for your example here. Yes, I can see in my data that some of the population-level posteriors are skewed -- typically towards the prior.
term "skew" has a technical meaning in statistics, and I have
a feeling you may be confusing it with the technical term "bias".
> This isn't true for all the cases where there is a difference however. See the attached for an example of relatively non-skewed posteriors -- the top is the population-level distribution, and the bottom are two of the individual study estimates. All the individual estimates are comparable to those that I show here, specifically none have their centres near the population-level estimates.
> This is an issue for me, because at the moment I can't really justify using the population-level estimates as a summary of the individual experiments. My experimental collaborators are -- quite rightly -- concerned that the population-level estimates aren't anything but a statistical mirage.
want to look at the posterior for the whole model. Assuming
that population level point estimates (you haven't said which
point estimates you're using---MLE, posterior median, posterior
mean) are representative of the population is also not right.
See, for example:
https://en.wikipedia.org/wiki/Ecological_fallacy
which in turn cites:
https://en.wikipedia.org/wiki/Simpson%27s_paradox
> The idea of my collaborators is just to use a non-hierarchical model, and just graph the individual estimates. I want to avoid this if possible, because I really see the value in this work as synthesising the data from previous works. However at the moment I think it might be the least bad option.
>
> Any thoughts?
then do the expedient thing for them and the editors for your career and
make it a longer-term goal to get the stats "right". I'm not just being
sarcastic here---it's what I have to do, too. See, e.g.,
https://lingpipe-blog.com/2014/10/29/beckys-and-my-annotation-paper-in-tacl/
- Bob
Bob Carpenter
Aug 25, 2016, 2:03:21 PM8/25/16
to stan-...@googlegroups.com
In case you didn't see, I went on to fill in a lot of direct
responses to points you brought up. Very good questions, by
the way. We really need to get our hands-on, applied stats book
written (and I promise no passives will be used, Michael :-)).
- Bob
responses to points you brought up. Very good questions, by
the way. We really need to get our hands-on, applied stats book
written (and I promise no passives will be used, Michael :-)).
- Bob
Ben Lambert
Aug 25, 2016, 2:15:37 PM8/25/16
to Stan users mailing list
On Thu, Aug 25, 2016 at 6:40 PM, Bob Carpenter <ca...@alias-i.com> wrote:
Yes, don't concentrate on point estimates.
Thanks for your patience here, Bob. And my apologies for lack of clarity. I am not concentrating on point estimates. I am looking at the whole posterior. Sometimes the bulk of the population-level probability mass is nowhere near the bulk of the mass for the individual estimates -- it is shifted. The argument remains the same regardless of whether I talk about a posterior summary interval -- like 10-90%, 25-75% quantiles -- or a point estimate. I am using intervals here to summarise my posteriors. The issue is that they do not necessarily overlap between the population-level estimates and the study-level ones.
I just can't shake the feeling that there is something odd going on in my code, or model. For all the other hierarchical models I've ever estimated the individual estimates are fairly-intuitively orientated in some way around the population-mean. Here this is the exception not the rule for each of the groups.
Also have answers inline, below.
- Bob
> On Aug 25, 2016, at 5:45 PM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:
>
> Hi Bob,
>
> Thanks for your example here. Yes, I can see in my data that some of the population-level posteriors are skewed -- typically towards the prior.
I still don't understand what you mean by "skewed" here. The
term "skew" has a technical meaning in statistics, and I have
a feeling you may be confusing it with the technical term "bias".
Sorry, yes I mean bias.
> This isn't true for all the cases where there is a difference however. See the attached for an example of relatively non-skewed posteriors -- the top is the population-level distribution, and the bottom are two of the individual study estimates. All the individual estimates are comparable to those that I show here, specifically none have their centres near the population-level estimates.
By "centre" you mean median?
I have used language -- mistakenly -- that gives too much emphasis on point estimates. I don't mean this. What I mean is that the posterior quantile-based intervals are not overlapping.
> This is an issue for me, because at the moment I can't really justify using the population-level estimates as a summary of the individual experiments. My experimental collaborators are -- quite rightly -- concerned that the population-level estimates aren't anything but a statistical mirage.
Treating point estimates with certainty is the mirage---you
want to look at the posterior for the whole model. Assuming
that population level point estimates (you haven't said which
point estimates you're using---MLE, posterior median, posterior
mean) are representative of the population is also not right.
See, for example:
https://en.wikipedia.org/wiki/Ecological_fallacy
Thanks for the reference. However I don't think this necessarily covers all of my cases. Particularly am thinking of the mean and median discussion in the reference. In many of them the mean, median, mode, posterior 10-90%/25-75% quantiles of my individual-study estimates are not near the population-level estimates.
which in turn cites:
https://en.wikipedia.org/wiki/Simpson%27s_paradox
> The idea of my collaborators is just to use a non-hierarchical model, and just graph the individual estimates. I want to avoid this if possible, because I really see the value in this work as synthesising the data from previous works. However at the moment I think it might be the least bad option.
>
> Any thoughts?
From a practical point of view, if you need to publish with collaborators,
then do the expedient thing for them and the editors for your career and
make it a longer-term goal to get the stats "right". I'm not just being
sarcastic here---it's what I have to do, too. See, e.g.,
https://lingpipe-blog.com/2014/10/29/beckys-and-my-annotation-paper-in-tacl/
Thanks Bob. No I take your point here.
- Bob
--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/itcTZmFAfFE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
Bob Carpenter
Aug 25, 2016, 3:22:21 PM8/25/16
to stan-...@googlegroups.com
Always a chance something's wrong with the code or the model
is misspecified. You can usually find problems in code if
you can't get the right coverage for your posterior intervals for
simulated data.
Could the problem be looking at the posteriors for the parameters
vs. the predictive posterior for new low-level parameters?
For example, if my hierarchical model for low-level parameters
beta[i] is
for (i in 1:I)
beta[i] ~ normal(mu, sigma);
then you don't necessarily expect the marginal posterior for mu
to overlap with the marginal posteriors for any of the beta[i].
What you want to look at is the posterior predictive distribution
for the lower-levle parameters, which in Stan looks like this:
generated quantities {
real beta_tilde;
beta_tilde = normal_rng(mu, sigma);
}
The posterior on beta_tilde is what you want to look at for
coverage, not the posterior over mu. The distribution for
beta_tilde draws will be
p(beta_new | y) = INTEGRAL_{mu, sigma} normal(beta_new | mu, sigma)
* p(mu, sigma | y)
d.mu d.sigma
where y here is all the original data. This distribution is
your prediction for the population based on the model and observed
data y.
- Bob
> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
is misspecified. You can usually find problems in code if
you can't get the right coverage for your posterior intervals for
simulated data.
Could the problem be looking at the posteriors for the parameters
vs. the predictive posterior for new low-level parameters?
For example, if my hierarchical model for low-level parameters
beta[i] is
for (i in 1:I)
beta[i] ~ normal(mu, sigma);
then you don't necessarily expect the marginal posterior for mu
to overlap with the marginal posteriors for any of the beta[i].
What you want to look at is the posterior predictive distribution
for the lower-levle parameters, which in Stan looks like this:
generated quantities {
real beta_tilde;
beta_tilde = normal_rng(mu, sigma);
}
The posterior on beta_tilde is what you want to look at for
coverage, not the posterior over mu. The distribution for
beta_tilde draws will be
p(beta_new | y) = INTEGRAL_{mu, sigma} normal(beta_new | mu, sigma)
* p(mu, sigma | y)
d.mu d.sigma
where y here is all the original data. This distribution is
your prediction for the population based on the model and observed
data y.
- Bob
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> For more options, visit https://groups.google.com/d/optout.
>
>
> 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 unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
Ben Lambert
Aug 25, 2016, 3:31:27 PM8/25/16
to Stan users mailing list
Thanks Bob. I am already looking at what you call 'beta_tilde' in my code. This is what differs from the individual study-level estimates. I do this using,
beta_tilde = normal_rng(a,sigma);
while(beta_tilde<0) beta_tile = normal_rng(a,sigma);
where I need the while loop to ensure that I get a positive value of this parameter.
The trouble is that the bulk of probability mass for a -- the population level parameter -- is sometimes nowhere near bulk of mass for the beta[i]s. This means that I get an estimated beta_tilde that is not close to the individual beta[i]s.
> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> 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+unsubscribe@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/itcTZmFAfFE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
Bob Carpenter
Aug 25, 2016, 5:12:34 PM8/25/16
to stan-...@googlegroups.com
Great. We need to get truncated RNGs---it's an issue that
nobody has time to work on. But the point is that it
doesn't matter if a isn't close if sigma is wide.
Normal priors can be pretty weak, so lots of individual-level
data will overwhelm the population estimates. If you
have a really spread out population, then you are still going
to get an estimate for a that may have a tight posterior that's
nowhere near the actual values. But in that case, sigma will
be large, so beta_tilde will be much wider and should still get
the right coverage.
Anyway, glad to see you're doing the right posterior predictive
checks here.
- Bob
> > To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
nobody has time to work on. But the point is that it
doesn't matter if a isn't close if sigma is wide.
Normal priors can be pretty weak, so lots of individual-level
data will overwhelm the population estimates. If you
have a really spread out population, then you are still going
to get an estimate for a that may have a tight posterior that's
nowhere near the actual values. But in that case, sigma will
be large, so beta_tilde will be much wider and should still get
the right coverage.
Anyway, glad to see you're doing the right posterior predictive
checks here.
- Bob
> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
> >
> >
> > --
> > 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.
> > For more options, visit https://groups.google.com/d/optout.
> >
> >
> > --
> > You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/itcTZmFAfFE/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/itcTZmFAfFE/unsubscribe.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> 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.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
Reply all
Reply to author
Forward
0 new messages