Varying-intercept & varying-slope model in rstan

[Maurits Evers]

Dear stan users.

I have been thoroughly enjoying working my way through McElreath’s and Gelman&Hill’s books. I have a question concerning the varying-intercept, varying-slope model, discussed in the radon example in Gelman&Hill (chapter 13), and the coffee robot example in McElreath (chapter 13).

I consider both individual level and group-level predictors (e.g. per-house radon and per-county uranium measurements), see attached stan model. 

McElreath and the Stan manual (section 5.12) provide a full hierarchical model with multivariate priors on the group-level coefficients. This involves a covariance matrix that accounts for the correlation between group-level intercepts and slopes, as given in my attached model. So far, so good.

1. Question 1: The rstan-transposed examples from Gelman&Hill give a much simpler model, with uncorrelated separate priors for the group-level intercepts and slopes (see the model here). Why the different modelling approaches? Intuitively I would assume that taking into account the correlation between intercepts and slopes allows for more efficient pooling of information across groups. So why the independent prior model in the stan examples for Gelman&Hill? 
2. Question 2: As an extension to the radon example, I would like to consider multiple uranium measurements, measured both at basement and ground floor levels of houses. I would like to explore whether some counties are showing unusually high or low radon levels, corrected for county-level uranium levels (and by taking into account variability in uranium levels between the basement and ground floor). Houses are not matched, i.e. we have uranium measurements of some houses within a county (basement and ground floor), and we have radon measurements of other houses within the same county (basement and ground floor). I am unsure how to model this within a hierarchical model.

Every time I feel I'm beginning to understand the big picture, I seem to be taking two steps back when faced with a new question. Any insight and explanation would be greatly appreciated!

Maurits

[Andrew Gelman]

Hi, it's better to allow corr between intercepts and slopes but in any case it's good to parameterize so that, a priori, you have no reason to expect a pos or neg correlation there.  My book with Jennifer was pre-Stan so we didn't have so many options for modeling and computation.

For q2:  The uranium measurements are not at the house level, they are at the county level.  It is radon measurements that are at the house level.  You can definitely have multiple radon measurements in a house; then you just need to include varying intercepts for houses, also you'll need coefs for whether the measurement is in basement or 1st floor.  We actually did this in our full data analysis which is discussed in our 1999 paper.

A

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

[Maurits Evers]

Dear Andrew.

Thanks for the fast reply. Could you please take a look at my inline comments below.

Thanks,

Maurits

On Tuesday, 8 November 2016 13:23:47 UTC+11, Andrew Gelman wrote:

Hi, it's better to allow corr between intercepts and slopes but in any case it's good to parameterize so that, a priori, you have no reason to expect a pos or neg correlation there.  My book with Jennifer was pre-Stan so we didn't have so many options for modeling and computation.

I understand. Just to clarify: I was referring to the stan-transposed examples from your book as they appear on github (https://github.com/stan-dev/example-models/wiki/ARM-Models). So are you saying that you did not model the correlation between slopes and intercepts to guarantee reproducibility with published lmer-based results? My confusion arose from Eq. (13.1) where in fact you do consider slope-intercept correlations.

 

For q2:  The uranium measurements are not at the house level, they are at the county level.  It is radon measurements that are at the house level.  You can definitely have multiple radon measurements in a house; then you just need to include varying intercepts for houses, also you'll need coefs for whether the measurement is in basement or 1st floor.  We actually did this in our full data analysis which is discussed in our 1999 paper.

Yes I understand. I was not interested in including multiple radon measurements per-house. Instead, my question was about modifying the model to include multiple uranium measurements per county (so let’s say from different houses again), both at basement and ground floor level. The question that I would like to address in my model would be: Is there a change in county-specific radon-levels between basement and ground floor, when we correct for basement and ground floor level uranium measurements? Do some counties exhibit an unusual change between basement and ground floor radon levels, independent of basement vs. ground floor changes in uranium levels. Does that make sense?

[Andrew Gelman]

Hi, I have not looked at the Stan-translated examples; they were done by someone else.  Jennifer and I are writing a new edition of our book (it will become 2 volumes) and I'll be doing all the Stan code then.  On your second question, no, it would not make sense to have multiple uranium measrurements per county.  Uranium measruements were taken from airplanes.  They are not at house level; they are just geographic measrurements that were then matched to the nearest county.  There is no such thing as a uranium measurement at the 1st floor or whatever.  Our 1999 paper may be helpful in explaining all this.

A

[Maurits Evers]

Hi Andrew.

Again, thanks for your quick response. I feel that we’re missing each other’s point. I am actually a nuclear physicist, so nuclear radiation measurements are something that I am very familiar with.

However, I may have been mis-using your radon example to explore a different model (and question). For the sake of the model’s theme, let’s assume that we have Polonium measurements (ignoring the uranium measurements entirely). Radon alpha-decays to polonium, and we have managed to measure the cumulative radiation of all polonium isotopes in different houses, at basement and ground floor levels, and in different counties.

My question is: How to model a floor-dependent effect of polonium radiation, whilst controlling for floor-dependent radon radiation; based on different per-house measurements, in different counties.

If I were to write down a model it would be something like this:

Here, j denote the county indices. 

I am not sure how to model this in rstan. The first two lines essentially correspond to two independent linear models. So one could estimate intercepts and slopes separately in two independent stan models, and then estimate coefficients c,d,e,f? In which case we wouldn't really propagate county- and individual measurement-based uncertainties across the two models within a Bayesian framework, or would we? 

Maurits

[Bob Carpenter]

If you write a joint model p(theta, phi, x, y) with

two sets of parameters theta and phi that factor as

  p(theta , phi | x, y) = p(theta | x, y) * p(phi | x, y)

then you'll get the same result fitting them together

as separately.  If they have overlapping parameters, on

the other hand, or one's conditional on the other, then the

info will flow both ways.  Your example does not look like

it'll factor because the a^polonium_j depends on a^radon_j.

I'd suggest generating fake data and fitting it to your

model to see if it beheaves as you expect.

- Bob

On Nov 7, 2016, at 11:15 PM, Maurits Evers <maurit...@gmail.com> wrote:

Hi Andrew.

Again, thanks for your quick response. I feel that we’re missing each other’s point. I am actually a nuclear physicist, so nuclear radiation measurements are something that I am very familiar with.

However, I may have been mis-using your radon example to explore a different model (and question). For the sake of the model’s theme, let’s assume that we have Polonium measurements (ignoring the uranium measurements entirely). Radon alpha-decays to polonium, and we have managed to measure the cumulative radiation of all polonium isotopes in different houses, at basement and ground floor levels, and in different counties.

My question is: How to model a floor-dependent effect of polonium radiation, whilst controlling for floor-dependent radon radiation; based on different per-house measurements, in different counties.

If I were to write down a model it would be something like this:

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+un...@googlegroups.com.

[Maurits Evers]

Hi Bob.

Yes, I agree in that the joint probability density does not factor into two independent probability densities. It that were the case we could fit them separately.

If I rewrite the code model as 

this becomes even clearer.  

The first line seems like a normal varying-intercept, varying-slope model with two additional group-dependent parameters that are determined from an independent model fit where we regress y^radon on x. 

I would like to double-check with you (and/or others), whether this model makes sense and — given adequate priors — could be fitted in stan as written here?

Best,

Maurits

PS. Yes, fake data will be my next step.

[Bob Carpenter]

Assuming the parameters make sense physically, it makes
sense statistically and should work the way you expect it to.
As long as you think generatively like this, building observed
and derived quantities bottom up from priors, you're in good
shape with Bayesian modeling. The priors can be uniform if
the posteriors are proper, but we recommend at least weakly informative
priors defining the expected order/scale of the answer.

Yes, you can code it in Stan directly as written with loops
or you can do it more efficiently with vectorization. Basically,
as long as you can write down a differentiable density (using
arithmetic and the functions we support), you can code it in Stan.

- Bob

> On Nov 9, 2016, at 7:49 PM, Maurits Evers <maurit...@gmail.com> wrote:
>
> Hi Bob.
>
> Yes, I agree in that the joint probability density does not factor into two independent probability densities. It that were the case we could fit them separately.
>
> If I rewrite the code model as
>
>

> this becomes even clearer.
>
> The first line seems like a normal varying-intercept, varying-slope model with two additional group-dependent parameters that are determined from an independent model fit where we regress y^radon on x.
>
> I would like to double-check with you (and/or others), whether this model makes sense and — given adequate priors — could be fitted in stan as written here?
>
> Best,
> Maurits
>
> PS. Yes, fake data will be my next step.
>
>
> On Thursday, 10 November 2016 03:21:53 UTC+11, Bob Carpenter wrote:
> If you write a joint model p(theta, phi, x, y) with
> two sets of parameters theta and phi that factor as
>
>
> p(theta , phi | x, y) = p(theta | x, y) * p(phi | x, y)
>
> then you'll get the same result fitting them together
> as separately. If they have overlapping parameters, on
> the other hand, or one's conditional on the other, then the
> info will flow both ways. Your example does not look like
> it'll factor because the a^polonium_j depends on a^radon_j.
>
>
>

[Maurits Evers]

Many thanks for your replies and advice, Bob (& Andrew).