Puzzlement over vectorization in BUGS epil example

[Larry Hunsicker]

Greetings, folks:  Back for more help.  I am having difficulty visualizing how the vectorized commands are working in the epil example (Book 1 of the old BUGS examples).  This is a mixed Poisson model.  (Ancient history, I know.  But the great masters tried to understand historical models before doing something "original.")

In this example, there are N (59) subjects, each seen T (4) times.  The dependent variable  y is declared in a 59 x 4 matrix.  Other covariates (constant within subject) are given in length 59 vectors.  There is an indicator variable for the fourth visit (all other visits are treated the same).  So visit number is given as a length four vector (0,0,0,1).  The data are eventually all centered, but that is, I think, irrelevant to my question.  b and b1 are two variance variables to be estimated for the patient and the patient * visit levels.  a0 and Alphas are all parameters to be fitted.  After various declarations and data manipulation, the final model is given as:

for(n in 1:N) {
    b[n] ~ normal(0, sigma_b);
    y[n] ~ poisson_log(a0 + alpha_Base * log_Base4_c[n]
                          + alpha_Trt * Trt_c[n]
                          + alpha_BT  * BT_c[n]   //log(Base) * trt
                          + alpha_Age * log_Age_c[n]
                          + b1[n] + alpha_V4 * V4_c + b[n]);
  }

My confusion comes in the last line.  There are 59 * 4 observations, with no missing values, and I am assuming that a (perhaps virtual) 59 * 4 matrix of predictors for each subject is being created.  I assume that a 59 * 4 length y vector is being created for each patient and each visit.  Since R and STAN go by column major order, the y vector holds the value for the first visit for each patient, then for the second visit for each patient, and so forth.  No problem with the vectors for the various columns of of patient level data.  Each column is just repeated four times, lining up correctly with the y's.  But the V4_c  variable is only length 4.  It should be repeated 59 times, but with as "rep(V4_c, each = 59), rather than as "rep(V4_c, 59) which would assign the indicators for the V4 visit incorrectly.  Since "V4_c" is referenced in the last line without any indices, how does the code know to create a column of "rep(v4_c, each = 59)" rather than rep(V4_c, 59), i.e., the length four V4c vector repeated 59 times? 

And am I visualizing correctly how the final design matrix and dependent variable vector are being created?

Thanks in advance for your help with my confusion.

Larry Hunsicker

[Larry Hunsicker]

OK.  I think that I understand it.  (I used the "print()" trick within 'epil.stan' to clarify what was meant by y(n).)   In STAN, one can refer to a whole row of y with y[n] (unlike R, where one would have to write "y[n,]").  So the four values of y[n] are each matched up with the covariate[n] scalar values, and also with, in turn, each of the four values of V4_t, which can be referenced without an index.   Is this correct?
Just out of curiosity, within the STAN code is a length N*T vector of y and a corresponding (N*T)*# of covariates +1 (the extra one for V4_t) design matrix ever created?
Again, thanks for your patience with me.
Larry Hunsicker

[Bob Carpenter]

> On Jun 16, 2015, at 10:39 AM, Larry Hunsicker <lhuns...@mchsi.com> wrote:
>
> OK. I think that I understand it. (I used the "print()" trick within 'epil.stan' to clarify what was meant by y(n).) In STAN, one can refer to a whole row of y with y[n] (unlike R, where one would have to write "y[n,]"). So the four values of y[n] are each matched up with the covariate[n] scalar values, and also with, in turn, each of the four values of V4_t, which can be referenced without an index. Is this correct?

Correct.

> Just out of curiosity, within the STAN code is a length N*T vector of y and a corresponding (N*T)*# of covariates +1 (the extra one for V4_t) design matrix ever created?

No, we don't need to create the whole design matrix in general.

The only requirement for Stan is that the model define the log density.

- Bob

[Larry Hunsicker]

Thanks, Bob.
Now one last question.  I used extract(model, pars = c('b1', 'b') with 4000 samples and confirmed that dim(b1) was (4000, 59) and dim(b) was (4000, 59, 4).  So, because b1 was declared outside of the for-loop, one slot for an estimate was created for each iteration of the for-loop (length(b1 = 59).  But because b was declared inside the for-loop, a slot was created for each estimate within each for-loop instance -- four, since each level of y had four values (and V4_c had four values) dim(b = (59,4)..  Do I have that right now?  If so, I think that I finally understand the code.

And thanks again for your support and rapid responses as I go through my learning curve.

Larry Hunsicker

[Bob Carpenter]

The quantities being estimated are fully determined by the declarations
in these blocks:

parameters
transformed parameters
generated quantities

The statements in the code, like loops, only define the density, not
the variables being estimated.

- 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.
> For more options, visit https://groups.google.com/d/optout.

[Larry Hunsicker]

Bob:
My apologies for being unclear in my last post.  Let me try to ask my question in a way that I hope is closer to correct.  b1 is intended to be the patient level residual and b the residual patient * visit residual, to account for extra Poisson variability.  Neither of these is an input value.  They are both declared to be draws from normal 0 mean distributions with sd's = sigma_b1 and sigma_b, respectively. 
My question is not how these values are fit by STAN, but only about how STAN knows that b1 is the patient level residual and b the patient by visit residual.  My guess is that this is because b1 is declared OUTSIDE of the do-loop, so that it sees that it has N levels, whereas b is declared INSIDE of the do-loop, so that it sees that, at each level of N, there are four (visit) levels.  I think that this is corroborated by the fact that at each draw, STAN draws N sample values for b1, but N * 4 values for b. 
At this point of my STAN learning curve, I don't need to know how STAN figures out the densities.  That comes later after I have made more progress through Gelman et al.  All I need to understand better today is how to code this sort of nested mixed model correctly.
Again, thanks for your help, your availability, and (I hope) your patience with me.
Larry

[Ben Goodrich]

On Tuesday, June 16, 2015 at 7:28:02 PM UTC-4, Larry Hunsicker wrote:

My question is not how these values are fit by STAN, but only about how STAN knows that b1 is the patient level residual and b the patient by visit residual.  My guess is that this is because b1 is declared OUTSIDE of the do-loop, so that it sees that it has N levels, whereas b is declared INSIDE of the do-loop, so that it sees that, at each level of N, there are four (visit) levels.  I think that this is corroborated by the fact that at each draw, STAN draws N sample values for b1, but N * 4 values for b.

What Bob was saying is that Stan (unlike BUGS) does not infer the number of elements in b and b1 from whether they are outside or inside some loop but simply accepts the dimensions that the user declares in the Stan program. In the epil example,

https://raw.githubusercontent.com/stan-dev/example-models/master/bugs_examples/vol1/epil/epil.stan

we have in the parameters block

...

real  b1[N]; 
vector[T] b[N];
...

so b1 has N elements and b has T * N elements regardless of what goes on in the subsequent blocks.

Ben

[Bob Carpenter]

The key to coding a model in a Stan program is to
understand how parameters are declared (with constraints)
and how the log probability function is defined (sampling
and increment_log_prob statements, plus implicit adjustment
for constraints). Everything else follows from that.

- Bob

[Larry Hunsicker]

AAAaaaaahhhhhh!!  Now I get it. (Or at least, I think that I have it now.)  I had either missed or not understood that declaration. b is a length T vector of N length vectors.  Then, since b1 is only N long, the draws for those numbers have to be the same for all 4 estimates of log(y) at each level of N (patient offset).  b has N  * T draws, so that the each of the four estimates of log(y) can have a separate draw (patient * visit offset).  . 

I had missed that, or perhaps misunderstood it.  "Now is the way clear!" (TS Eliot)

>> The key to coding a model in a Stan program is to understand ... how the log probability function is defined (sampling and increment_log_prob statements, plus implicit adjustment  for constraints).

That part is for next week... or next year, or something.

Thanks to the two of you for your patience and persistence.

Larry

[Larry Hunsicker]

Bob and Ben (or anyone else that can help me):

OK.  I apologize.  But I have another question. 

Based on yesterday's discussion, I thought that I understood that the prior for b: "b[n] ~ normal(0, sigma_b);" could be put pretty much anywhere, since the shape of b had been set by its earlier declaration in the parameters section.  But when I moved the prior for b out of the do-loop, the compile failed with a message to the effect that the type parameter for the object didn't match up with the earlier declaration.  So I looked up the parameter types and tried
"b[N][K] ~ normal(0, sigma_b)".  This compiled and ran, but glacially slowly.  I lost patience and cancelled the call, so I don't know whether I would have gotten the same answer as when the prior for b is declared inside of the do-loop.  But either that change was an out and out mistake, or the placement of the prior led to an awful de-optimization.  So it is clear that there is still something that I don't understand correctly about why the prior was placed inside of the do-loop.

Why is it important to put the prior for b[n] inside of the do-loop?  I could, of course, just use epil.stan as a template for future analyses of this sort.  But I'd much prefer to understand what I am doing.  Again, if this is all explained somewhere in Gelman et. al, I'd be happy just to receive a pointer to its location.

Sorry to be such a pest.  But you two have been so helpful...

Larry

[Ben Goodrich]

On Wednesday, June 17, 2015 at 6:56:10 PM UTC-4, Larry Hunsicker wrote:

Based on yesterday's discussion, I thought that I understood that the prior for b: "b[n] ~ normal(0, sigma_b);" could be put pretty much anywhere, since the shape of b had been set by its earlier declaration in the parameters section. 

It is true that the dimensions of b are declared in the parameters section, but a for-loop in Stan

  for(n in 1:N) {

    b[n] ~ normal(0, sigma_b);

 
    ...

works like a for-loop in any other language. Here it indicates that the n-th element of b --- which happens to be a vector of length T in this case --- has a prior on its elements that is iid normal with mean zero and standard deviation sigma_b.

• If you put anything involving b[n] outside the for loop, Stan wouldn't know what n refers to.
• b[N][K] ~ normal(0, sigma_b) does not make much sense because it puts the normal(0, sigma_b) prior only on the N,K element and does not put a proper prior on any of the other elements (which is why it took a long time to sample)

There are some ways that this could be improved upon:

• Delete the line that says b[n] ~ normal(0, sigma_b); and add (outside the for loop) the line to_array_1d(b) ~ normal(0, sigma_b);. This concatenates all the elements of b into a one-dimensional array and asserts the prior is iid normal with mean 0 and standard deviation sigma_b.
• Delete the line that says b[n] ~ normal(0, sigma_b); and the line that says b1 ~ normal(0, sigma_b1); Add (outside the for loop) the line to_array_1d(b) ~ normal(0, 1); and the line b1 ~ normal(0,1). Finally, change the last line of the likelihood to say + b1[n] * sigma_b1 + alpha_V4 * V4_c + b[n] * sigma_b); . This basically redefines b[n] and b1 to be standard normal shifts that are scaled by the unknowns sigma_b and sigma_b1 respectively, which often improves the efficiency of the sampling by reducing the dependence between these unknown parameters.

You can always print things to the screen if you are having trouble understanding the structure:

print("b is");
for(n in 1:N) print(b[n]);

but something like that is best called in diagnostic mode (in R with test_grad = TRUE and chains = 1) so that you only see it once rather than millions of times.

Ben

[Larry Hunsicker]

Many thanks, Ben.  That is very clear.  I'm going to try out your alternatives and time them against the original.  But you have really helped me understand this business.

>>  You can always print things to the screen if you are having trouble understanding the structure:

Yes.  I have tried this at the end of the transformed variables section to be sure that I had got the arithmetic correct.

>>  but something like that is best called in diagnostic mode (in R with test_grad = TRUE and chains = 1) so that you only see it once rather than millions of times.

And yes.  I discovered this, too.  I had to kill my R instance to stop the torrent of red.  I'll try your suggestion above.

I think that that is, finally, the end of this thread.  Many, many thanks.

Larry Hunsicker

[Larry Hunsicker]

Ben:
Back again, because I can't make your first suggested revision work.  The definition of b from the parameters block is:
  vector[T] b[N];
The model definition from my version of the epil model statement (that works) is:

model {
  a0 ~ normal(0, 100);
  alpha_Base ~ normal(0, 100);
  alpha_Trt  ~ normal(0, 100);
  alpha_BT   ~ normal(0, 100);
  alpha_Age  ~ normal(0, 100);
  alpha_V4   ~ normal(0, 100);
  sigmasq_b1 ~ inv_gamma(.001, .001);
  sigmasq_b ~ inv_gamma(.001, .001);
  b1 ~ normal(0, sigma_b1);

  for(n in 1:N) {
    b[n] ~ normal(0, sigma_b);

    y[n] ~ poisson_log(a0 + alpha_Base * lBase4c[n]
                       + alpha_Trt * trtc[n]
                       + alpha_BT  * BTc1[n]
                       + alpha_Age * lagec[n]
                       + alpha_V4 * V4c
                       + b1[n] + b[n]);
  }
}

I followed the first of your suggestions; the only the changes were to the bottom half of the model:
...
  sigmasq_b ~ inv_gamma(.001, .001);
  b1 ~ normal(0, sigma_b1);
  to_array_1d(b) ~ normal(sigma_b);
  for(n in 1:N) {
//    b[n] ~ normal(0, sigma_b);
    y[n] ~ poisson_log(a0 + alpha_Base * lBase4c[n]
                       + alpha_Trt * trtc[n]
                       + alpha_BT  * BTc1[n]
...

This failed to compile with the following error report:

Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
  failed to parse Stan model 'epil3' with error message:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
no matches for function name="to_array_1d"
    arg 0 type=vector[1]
available function signatures for to_array_1d:
0.  to_array_1d(matrix) : real[1]
1.  to_array_1d(vector) : real[1]
2.  to_array_1d(row vector) : real[1]
3.  to_array_1d(real[1]) : real[1]
4.  to_array_1d(int[1]) : int[1]
5.  to_array_1d(real[2]) : real[1]
6.  to_array_1d(int[2]) : int[1]
7.  to_array_1d(real[3]) : real[1]
8.  to_array_1d(int[3]) : int[1]
9.  to_array_1d(real[4]) : real[1]
10.  to_array_1d(int[4]) : int[1]
11.  to_array_1d(real[5]) : real[1]
12.  to_array_1d(int[5]) : int[1]
13.  to_array_1d(real[6]) : real[1]
14.  to_array_1d(int[6]) : int[1]
15.  to_array_1d(real[7]) : real[1]
16.  to_array_1d(int[7]) : int[1]
17.  to_array_1d(real[8]) : real[1]
18.  to_array_1d(int[8]) : int[1]
19.  to_array_1d(real[9]) : real[1]
20.  to_array_1d(int[9]) : int[1]
no matches for function name="to_array_1d"
    arg 0 type=vector[1]
available function signatures for to_array_1d:
0.  to_array_1d(matrix) : real[1]
1.  to_array_1d(vector) : real[1]
2.  to_array_1d(row vector) : real[1]
3.  to_array_1d(real[1]) : real[1]
4.  to_array_1d(int[1]) : int[1]
5.  to_array_1d(real[2]) : real[1]
6.  to_array_1d(int[2]) : int[1]
7.  to_array_1d(real[3]) : real[1]
8.  to_array_1d(int[3]) : int[1]
9.  to_array_1d(real[4]) : real[1]
10.  to_array_1d(int[4]) : int[1]
11.  to_array_1d(real[5]) : real[1]
12.  to_array_1d(int[5]) : int[1]
13.  to_array_1d(real[6]) : real[1]
14.  to_array_1d(int[6]) : int[1]
15.  to_array_1d(real[7]) : real[1]
16.  to_array_1d(int[7]) : int[1]
17.  to_array_1d(real[8]) : real[1]
18.  to_array_1d(int[8]) : int[1]
19.  to_array_1d(real[9]) : real[1]
20.  to_array_1d(int[9]) : int[1]
expression is ill formed

ERROR at line 53
 51:      sigmasq_b ~ inv_gamma(.001, .001);
 52:      b1 ~ normal(0, sigma_b1);
 53:      to_array_1d(b) ~ normal(sigma_b);
                         ^
 54:      for(n in

I think that this is tellling me that the function "to_array_1d" doesn't have a version that knows how to handle a vector of real vectors.  I fished around, trying to_array_2d, and tried declaring b to be a matrix or a two dimensional array.  But none of this worked. 

Any suggestions about this? 

Your second suggestion, to redefine b[n] and b1, compiles, runs, and gives a correct answer. 

b1 ~ normal(0, 1);
  for(n in 1:N) {
    b[n] ~ normal(0, 1);
    y[n] ~ poisson_log(a0 + alpha_Base * lBase4c[n]
                       + alpha_Trt * trtc[n]
                       + alpha_BT  * BTc1[n]
                       + alpha_Age * lagec[n]
                       + alpha_V4 * V4c
                       + b1[n] * sigma_b1 + b[n] * sigma_b);
  }

But it is actually much slower than the original code:
Old:
   user  system elapsed
  83.88    0.02  115.11
New:
  user  system elapsed
 108.11    0.22  145.84

I have no clue about why this should be.

Again, thanks with your patience with me. 

Larry Hunsicker

[Andrew Gelman]

Umm, I can’t keep doing this, but:

- Don’t ever do inv_gamma(.001,.001)

- Please scale your parameters to be roughly on unit scale, in which case don’t use normal(0,100) as a noninformative prior.  normal(0,10) or even normal(0,1) will do the trick.

Also the model below seems wack because b1 and b do not seem separately identified.

[Larry Hunsicker]

Thanks, Andrew:

The code is all based on the old BUGS examples, and all the BUGS examples seem to use those very wide priors.  I have fiddled with the normal priors and agree that normal(0, 10), or even normal(0,2) seem to give very similar posteriors to normal(0, 100) for, e.g. sd of average height of women in meters.  But I haven't got to learning about sensible non- or weakly- informative priors yet.  Maybe I should just switch over to the ARM examples to learn from.  I promise not to use inv_gamma(.001, .001) ever again. 

Larry Hunsicker

[Andrew Gelman]

Thx.  And inv_gamma(.001,.001) is not wide at all!

[Bob Carpenter]

Sorry -- I don't have time to carefully go through all the
syntax of everyone's model.

> Any suggestions about this?

Use a matrix. And to_vector. Or just leave it as is.

> Your second suggestion, to redefine b[n] and b1, compiles, runs, and gives a correct answer.
>
> b1 ~ normal(0, 1);
> for(n in 1:N) {
> b[n] ~ normal(0, 1);
> y[n] ~ poisson_log(a0 + alpha_Base * lBase4c[n]
> + alpha_Trt * trtc[n]
> + alpha_BT * BTc1[n]
> + alpha_Age * lagec[n]
> + alpha_V4 * V4c
> + b1[n] * sigma_b1 + b[n] * sigma_b);
> }
>
> But it is actually much slower than the original code:
> Old:
> user system elapsed
> 83.88 0.02 115.11
> New:
> user system elapsed
> 108.11 0.22 145.84
>
> I have no clue about why this should be.

What's b1? You don't want to put b inside and outside the loop.

There's a lot of variation between runs, so I wouldn't
count that as huge --- it might just be noise.

- Bob