Stan code for Gelman and Hill (2003)

[Thomas Zumbrunn]

Dear (R)Stan users

I was wondering if anyone has translated or knows of someone who has
translated all the BUGS code in

Gelman A & Hill J (2006) Data Analysis Using Regression and
Multilevel/Hierarchical Models. Cambridge University Press

to Stan code. I'm reading the book with a couple of colleagues, and we'd like
to use Stan instead of (Open|Win)BUGS. Any hints are appreciated.

Thanks,
/thomas

[Tim Triche, Jr.]

Not unrelated, will Bayesian Data Analysis 3ed be written with examples in Stan?  Because that would be awesome.

--
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/groups/opt_out.

--
A model is a lie that helps you see the truth.

Howard Skipper

[Andrew Gelman]

Hi, yes, BDA3 has new sections on HMC and Stan, and the revised Appendix C features Stan code.  We don't have Stan implenetations for all the BDA examples, but I wouldn't stop anyone from doing that themselves…

Also, we have not yet translated the Bugs code in ARM into Stan, but we plan to do so for the next edition.

A

[Bob Carpenter]

The new appendix is about Stan. But there's also
lots more in the text, like variational inference,
HMC, and non-parametric modeling.

I don't know that anyone's systematically translated
the examples for Gelman and Hill. We'd be delighted
to help someone do that if Andrew and Jennifer could find
the data for us. It's been on the to-do list from the
get go.

So if there are particular models you'd like help translating,
feel free to send them to the list and we'd be glad to
help out.

We'd also be happy to distribute the results with Stan.

- Bob

On 4/4/13 2:43 PM, Tim Triche, Jr. wrote:
> Not unrelated, will Bayesian Data Analysis 3ed be written with examples in Stan? Because that would be awesome.
>
>
>
> On Thu, Apr 4, 2013 at 11:39 AM, Thomas Zumbrunn <t.zum...@unibas.ch <mailto:t.zum...@unibas.ch>> wrote:
>
> Dear (R)Stan users
>
> I was wondering if anyone has translated or knows of someone who has
> translated all the BUGS code in
>
> Gelman A & Hill J (2006) Data Analysis Using Regression and
> Multilevel/Hierarchical Models. Cambridge University Press
>
> to Stan code. I'm reading the book with a couple of colleagues, and we'd like
> to use Stan instead of (Open|Win)BUGS. Any hints are appreciated.
>
> Thanks,
> /thomas
>
> --
> 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 <mailto:stan-users%2Bunsu...@googlegroups.com>.

> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
>
>
> --

> /A model is a lie that helps you see the truth./
> /
> /
> Howard Skipper <http://cancerres.aacrjournals.org/content/31/9/1173.full.pdf>

[Thomas Zumbrunn]

Thanks. I'm looking forward to getting my hands on BDA3.

/thomas

[Thomas Zumbrunn]

Thanks a lot. I'll give it a try and ask on the list if I run into problems
(which I will for sure).

/thomas



PS: Just realised that my original subject line was wrong - the year should've
been 2006 instead of 2003...

[Bob Carpenter]

Section 11.6 of the Stan manual has an example of
the basic ordered logistic model.

I'll have to refresh my memory on the storable votes
model. It shouldn't be too hard, but it might be
a couple of days until I can get to it. So please
remind me if I don't reply in a couple of days.

- Bob

On 4/9/13 11:23 PM, Duncan wrote:
> Well, since you've so kindly offered, Bob, I'd very much appreciate seeing the model under "17.6 Multilevel ordered
> categorical regression" (p. 383 of Gelman and Hill) translated into Stan!
> Duncan

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

[sp_r...@yahoo.it]

I've just started studying Bayesian methods, so I'm translating BUGS _and_ regression code in Gelman & Hill to Stan code.

I attach a file containing code from chapters 3 and 4.

I hope that my code is not too ugly, and that you could give me some advice. For example: is there a standard way to code discrete predictors (factors)? I've found my way when translating Ch.4/4.5_Other transformations.R, the output is consistent, but I'm not sure that my way is a clever one... This is why I didn't try to translate Ch.4/4.7_Fitting a series of regressions.R.

Sergio

[Bob Carpenter]

Thanks for sharing. We have an effort going there to
implement all of the Gelman and Hill models that those actually
doing it can fill you in on (Alexeij? Daniel?).

We plan to use our optimizers to recreate the functionality of
lmer (bglmer) max marginal estimates in Stan. And of course
to get all the BUGS models working.

I'd highly recommend separating three things by file:

1. Stan model

2. data set

3. R code for fitting model & plotting results

By having 1 and 2 in separate files, Stan users who don't use R
will benefit. And it makes it easier to browse the model
implementations.

As is, you have multiple models in a single R file.

Your code should run out of the box, either by including
copies of the data or using URLs. As is, you have the
URL code commented out and a dangling local reference.

In general, you shouldn't have commented out alternative code
included with the actual code --- it's just confusing for users.
If there are two ways to do thing, put them both in the code
and have some way to configure them. If not, remove the one
that's not needed.

There's also the choice of whether to do data transformations in
Stan or in the calling program (R in this case). Either the transform
should be done to create the data set that's included or it should
be done in Stan in order to preserve the standalone nature of
model+data.

- Bob

[sp_r...@yahoo.it]

Something like this?

Sergio

EDIT: I've removed my attachment. Please, see my next message.

[sp_r...@yahoo.it]

I've read a lot of messages and noticed that:
a) you use stopifnot(require(rstan)) instead of library(rstan);
b) another developer (I can't remember his name...) appends .RData to saved compiled models.
Moreover,
c) as far I can understand, when using Stan from the command line the data file can contain more variables than the ones declared in the data block, so I enclose just one data file, even if only a subset of its variables are meant to be read by a model;
d) remembering your words (my R code should run out of the box) I've rewritten my scripts; I'm used to execute a script line by line (or region by region) in Emacs+ESS, but you can now type source(<script name>) and look at the numerical and graphical results.

So I've removed my previous attachment (GH2007_Ch.3.tar.gz) and attach a new one: GH2007_Ch.3_0.2.tar.gz

Best wishes,

Sergio

[Bob Carpenter]

On 7/14/13 12:30 PM, sp_r...@yahoo.it wrote:
> I've read a lot of messages and noticed that:
> a) you use stopifnot(require(rstan)) instead of library(rstan);

I'll leave that to the R gurus. I use library(rstan).

> b) another developer (I can't remember his name...) appends .RData to saved compiled models.

The preferred form is .data.R if the data is in the
form of R source code (like a dump) or .csv if it's
a CSV file.

> Moreover,
> c) as far I can understand, when using Stan from the command line the data file can contain more variables than the ones
> declared in the data block, so I enclose just one data file, even if only a subset of its variables are meant to be read
> by a model;

Correct.

> d) remembering your words (my R code should run out of the box) I've rewritten my scripts; I'm used to execute a script
> line by line (or region by region) in Emacs+ESS, but you can now type source(<script name>) and look at the numerical
> and graphical results.

That's much easier for others to replicate.

> So I've removed my previous attachment (GH2007_Ch.3.tar.gz) and attach a new one: GH2007_Ch.3_0.2.tar.gz

Thanks.

- Bob

[Sergio Polini]

>> b) another developer (I can't remember his name...) appends .RData to
>> saved compiled models.
>
> The preferred form is .data.R if the data is in the
> form of R source code (like a dump) or .csv if it's
> a CSV file.

Sorry, my english is poor, but I was speaking of compiled models, not
data files. I.e:

sm <- stan_model(...)
save(sm, file="sm.RData")

then, when you need the compiled model again:

load("sm.RData", verbose=TRUE)

As to the data file (only one in the attached .tar.gz) its name is
kidhsiq.data.R.

Sergio

[Bob Carpenter]

On 7/14/13 1:53 PM, Sergio Polini wrote:
>>> b) another developer (I can't remember his name...) appends .RData to
>>> saved compiled models.
>>
>> The preferred form is .data.R if the data is in the
>> form of R source code (like a dump) or .csv if it's
>> a CSV file.
>
> Sorry, my english is poor, but I was speaking of compiled models, not data files. I.e:

Your English is fine --- I just didn't read closely enough, because
I wasn't even aware you can do this. I usually play with my models
enough that saving them is rarely an issue.

For R binaries, I think .Rdata is the right suffix.

I have no idea how stable the saved version of models is in R.
Will they break in the next R release? No big deal even if they
do because you can always recompile them if you keep the source
Stan model.

- Bob

[Peter Li]

The code you've attached looks great. I'm actually working on writing
the ARM models in Stan too. What I've done so far is in
https://github.com/stan-dev/stan/tree/feature/ARM. I need to do a bit of
reorganizing, but I'll mostly be following a similar format/organization
as you.

Alexij is also working on the ARM models, and from what I understand,
he's already done the majority of ~ch16-20 and I've got a fair amount of
the first 12 chapters done (without the R code). In any case, if you
want to work together on this let me know and we can coordinate who's
doing what so we don't end up writing up the same models several times. :-)


Peter

On 7/14/13 12:30 PM, sp_r...@yahoo.it wrote:

[sp_r...@yahoo.it]

I'ld say that my code looks decent since I've followed Bob's advice, and that can be improved. E.g., I'ld now use read_rdump instead of source.
I'm a learner. I'm learning Bayesian statistics, I'm learning (R)Stan (and I'm writing my bachelor's thesis ... about multilevel models in longitudinal/panel data analysis.)

So let me know what you think I could do, and I'll try to do it.

Sergio

[Peter Li]

On 7/16/13 6:17 PM, sp_r...@yahoo.it wrote:

I'ld say that my code looks decent since I've followed Bob's advice, and that can be improved. E.g., I'ld now use read_rdump instead of source.
I'm a learner. I'm learning Bayesian statistics, I'm learning (R)Stan (and I'm writing my bachelor's thesis ... about multilevel models in longitudinal/panel data analysis.)

That's a lot to juggle! I am also learning how to actually use Stan and about Bayesian inference. Coding these has been a great learning experience.

So let me know what you think I could do, and I'll try to do it.

You can work on whatever model/chapters you want really. I've got most of the stan code and data files for most of the models in the first 12 chapters of ARM in https://github.com/stan-dev/stan/tree/feature/ARM/src/models/ARM/models%20coded%20but%20need%20to%20organize, but I haven't done the R code yet, so you could do that if you want. Overall, it's really up to you. Just let me know, and I'll work on whatever you're not doing.

Also, Bob wants us to use ggplot2 instead of base graphics (see https://github.com/stan-dev/stan/tree/feature/ARM/src/models/ARM/Ch.3).


Peter

[Bob Carpenter]

On 7/16/13 6:42 PM, Peter Li wrote:
>
> ...

> Also, Bob wants us to use ggplot2 instead of base graphics (see
> https://github.com/stan-dev/stan/tree/feature/ARM/src/models/ARM/Ch.3).

That wasn't meant as a demand! It would be my preference, but
you all are driving this bus.

I do suggest you ask Andrew what his preference would be.

I think he and Jennifer will use all of this for the next edition
of their book, so it's more important to make Andrew happy than
me on this one!

- Bob

P.S. I also learned applied stats going through Gelman and Hill.
It's a great source for learning modeling. It also revs one up
to the point where Bayesian Data Analysis makes sense. :-)

[Andrew Gelman]

I have no strong feelings on graphs. I think that ggplot2 is preferred by the experts. We used base graphics in the rstan output because ggplot2 can be difficult to customize, and for the rstan output we had to cram lots of graphics into a small space. For the arm graphs, just use whatever works best for you.
A

[Peter Li]

I'll go with ggplot2 then. It's worked pretty well so far--just finished
adding code for ch3.


Peter

[Jiqiang Guo]

The ggplot2 ebook is free, or free throught Columbia library.  

http://link.springer.com/book/10.1007/978-0-387-98141-3/page/1

As Andrew mentioned, RStan does not use ggplot2.

--

Jiqiang 

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

For more options, visit https://groups.google.com/groups/opt_out.

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

[Sergio Polini]

Il 17/07/2013 00:42, Peter Li ha scritto:
> Overall, it's really up to you. Just let me know, and I'll work on
> whatever you're not doing.

It's very kind of you, but I think that a few guidelines would be useful.
I've seen that you have modified my code. If I write the code as you
like (not an issue, really) you save time. It's a matter of productivity ;-)

Yes, I can go on and write the R code. My first guidelines (Bob's
advice) are:
1) model and (stan_rdumped) data in separate files;
2) R code should run out of the box (source the script and look at
numerical and graphical results);
3) no commented out alternative code;
4) data transformation in Stan (transformed data block).

Other choices:
5) library() or stopifnot(require())?
6) saving and loading a compiled model;
7) source() or read_rdump()?
8) graphics.

My personal opinions:

5) stopifnot() could be used when a package is not part of a standard R
istallation (such as rstan, but also ggplot2!), but it could be
confusing; I'ld now leave sopifnot() to the R gurus and use library()
(Bob's advice);

6) I've used "if (!file.exists("sm.RData")) then compile&save else load"
because it's convenient, but I got a few segfaults when loading and
reloading several times the same DSO (see:
https://groups.google.com/d/msg/stan-users/0Ev3agQB4g0/pCYN2bwyabQJ).
So I'm going tu use "if (!exists(sm)) then { if
(file.exists("sm.RData")) then load else compile&save }".

7) read_rdump() is very nice, because one can write:

fit <- sampling(sm, data=read_rdump("foo.data.R"))

instead of

source("foo.data.R")
dataList <- list(...)
fit <- sampling(sm, data=dataList)

but looking at the variables could be a bit cumbersome because
read_rdump() does not read the data into the user's workspace.

8) The question is: who is the target user?
No problem if the answer is "whoever".
But if the answer was "whoever is reading ARM" then the R code should
just show how to use stan() instead of lm()/lmer()/bugs(), extract()
instead of sim(). BTW, sim() is part of the arm package, which loads
lattice :-)
Having to "think in ggplot2" instead of "think in base graphics" could
overload the target user.

However these are just my humble opinions.
If you favor, say, "sopifnot, shorter if, read_rdump, ggplot2", I'll
adapt to save your time.

Sergio

[Jiqiang Guo]

If function source is used, variable names can be specified instead of a list.  For example, 

source('foo.data.R')

stan(file = 'foo.stan', data = c('x', 'y')) 

--

Jiqiang 

[Tim Triche, Jr.]

If you use require() instead of library, you will get most of the functionality of stopifnot(library(foo)), and a user missing the library in question will be informed the lack of the library is the source of their problems.

IMHO this is the way to go, especially until rstan ends up on CRAN.

Sergio

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

For more options, visit https://groups.google.com/groups/opt_out.

--

A model is a lie that helps you see the truth.

Howard Skipper

[Peter Li]

On 7/17/13 6:24 AM, Sergio Polini wrote:
> Il 17/07/2013 00:42, Peter Li ha scritto:
>> Overall, it's really up to you. Just let me know, and I'll work on
>> whatever you're not doing.
>
> It's very kind of you, but I think that a few guidelines would be useful.
> I've seen that you have modified my code. If I write the code as you
> like (not an issue, really) you save time. It's a matter of
> productivity ;-)
>
> Yes, I can go on and write the R code. My first guidelines (Bob's
> advice) are:
> 1) model and (stan_rdumped) data in separate files;
> 2) R code should run out of the box (source the script and look at
> numerical and graphical results);
> 3) no commented out alternative code;
> 4) data transformation in Stan (transformed data block).

These all sound great. I may need to modify some of what I've done
already, but that's fine.

>
> Other choices:
> 5) library() or stopifnot(require())?
> 6) saving and loading a compiled model;
> 7) source() or read_rdump()?
> 8) graphics.
>
> My personal opinions:
>
> 5) stopifnot() could be used when a package is not part of a standard
> R istallation (such as rstan, but also ggplot2!), but it could be
> confusing; I'ld now leave sopifnot() to the R gurus and use library()
> (Bob's advice);

I have no opinion, so we can go with library if that's what you think is
best. I'm still rather new to R, so I'm learning a lot as I go. :)

>
> 6) I've used "if (!file.exists("sm.RData")) then compile&save else
> load" because it's convenient, but I got a few segfaults when loading
> and reloading several times the same DSO (see:
> https://groups.google.com/d/msg/stan-users/0Ev3agQB4g0/pCYN2bwyabQJ).
> So I'm going tu use "if (!exists(sm)) then { if
> (file.exists("sm.RData")) then load else compile&save }".

I'll need to modify my code to follow this. I haven't run into that
error yet, but I've run into another error where the code doesn't
compile. It says that I'm mixing types, but I haven't looked into it. It
shouldn't be too big of a problem.

>
> 7) read_rdump() is very nice, because one can write:
>
> fit <- sampling(sm, data=read_rdump("foo.data.R"))
>
> instead of
>
> source("foo.data.R")
> dataList <- list(...)
> fit <- sampling(sm, data=dataList)
>
> but looking at the variables could be a bit cumbersome because
> read_rdump() does not read the data into the user's workspace.

I sorta like the latter method (not read_rdump) because sometimes you
want to follow the same model, but transform the params a bit (i.e.
centering or something), so this way you don't have to write a whole new
stan file. But I guess this goes back to using the data transformation
block. Regardless, I've coded up the transformations in separate stan
files but do not call those in a lot of the R code--I can change this
easily though.

>
> 8) The question is: who is the target user?
> No problem if the answer is "whoever".
> But if the answer was "whoever is reading ARM" then the R code should
> just show how to use stan() instead of lm()/lmer()/bugs(), extract()
> instead of sim(). BTW, sim() is part of the arm package, which loads
> lattice :-)
> Having to "think in ggplot2" instead of "think in base graphics" could
> overload the target user.

I think the target is anyone reading ARM. I agree that we should be
using stan instead of lm/lmer/bugs. But I'm not sure about sim() and how
it works (same goes for extract()) so I'm not sure how to replace sim()
with extract(). If you could write up an example, that'd be great; then
I could just follow that.

I'm also unsure about the graphics packages, but supposedly ggplot2 is
prefered by the "experts", so I was just going to go with that. I mean,
the users probably aren't going to be reading the book to follow how the
graphics are coded. Also, if you're not comfortable with using ggplot2,
I can take care of that.


Peter

[Bob Carpenter]

On 7/17/13 2:44 PM, Peter Li wrote:
>
> On 7/17/13 6:24 AM, Sergio Polini wrote:
>> Il 17/07/2013 00:42, Peter Li ha scritto:

>> 8) The question is: who is the target user?
>> No problem if the answer is "whoever".
>> But if the answer was "whoever is reading ARM" then the R code should just show how to use stan() instead of
>> lm()/lmer()/bugs(), extract() instead of sim(). BTW, sim() is part of the arm package, which loads lattice :-)
>> Having to "think in ggplot2" instead of "think in base graphics" could overload the target user.

> I think the target is anyone reading ARM. I agree that we should be using stan instead of lm/lmer/bugs. But I'm not sure
> about sim() and how it works (same goes for extract()) so I'm not sure how to replace sim() with extract(). If you could
> write up an example, that'd be great; then I could just follow that.

We have two target audiences:

1. Stan users who want example models.

2. Readers of ARM who want to see the code.

We don't want to depend on the previous ARM package on CRAN.
We won't be able to put this on CRAN because of the dependency
on Stan (which is too big and too slow to build).

You definitely want to be using extract() in RStan to get
the samples out.

If sim() has to do with forward simulation (just a wild guess),
then you should do that in the generated quantities block of
Stan. That's why we include the RNGs that Peter wrote!

- Bob

[Sergio Polini]

Il 17/07/2013 20:44, Peter Li ha scritto:
> I sorta like the latter method (not read_rdump) because sometimes you
> want to follow the same model, but transform the params a bit (i.e.
> centering or something), so this way you don't have to write a whole new
> stan file. But I guess this goes back to using the data transformation
> block. Regardless, I've coded up the transformations in separate stan
> files but do not call those in a lot of the R code--I can change this
> easily though.

If "Stan users who want example models" are target audience, perhaps the
data should be transformed using the tranformed data block.

> I think the target is anyone reading ARM. I agree that we should be
> using stan instead of lm/lmer/bugs. But I'm not sure about sim() and how
> it works (same goes for extract()) so I'm not sure how to replace sim()
> with extract(). If you could write up an example, that'd be great; then
> I could just follow that.

Caveat: I'm a learner :-) However...

sim() gets posterior simulations of sigma and beta from a 'lm', of beta
from a 'glm' or 'mer' object.
If you look at the surce code, you see that sim() uses frequentist point
estimates, their covariance and R functions such as mvrnorm() and
rchisq() to simulate a posterior distribution.

Stan's output _is_ a (sampled) posterior distribution. If you use
uninformative priors and normal likelihood functions, you get similar
results.

I attach an example, ARM Figure 4.1 right.

> I'm also unsure about the graphics packages, but supposedly ggplot2 is
> prefered by the "experts", so I was just going to go with that. I mean,
> the users probably aren't going to be reading the book to follow how the
> graphics are coded. Also, if you're not comfortable with using ggplot2,
> I can take care of that.

ggplot2 is fine.

Sergio

[Bob Carpenter]

Thanks for the attachment.

Please keep the models separate. I don't like them in the
R code. We'll have Stan users using the command-line and our
Python interface, and putting the models into R makes it very
hard on them.

The following should all be in separate files:

- model

- data set

- R code

I know that's not how the first RStan example goes on the RStan
install guide, but that's just to keep things simple, not for
code you'd like non-R-users to also use.

On 7/17/13 6:13 PM, Sergio Polini wrote:
> Il 17/07/2013 20:44, Peter Li ha scritto:
>> I sorta like the latter method (not read_rdump) because sometimes you
>> want to follow the same model, but transform the params a bit (i.e.
>> centering or something), so this way you don't have to write a whole new
>> stan file. But I guess this goes back to using the data transformation
>> block. Regardless, I've coded up the transformations in separate stan
>> files but do not call those in a lot of the R code--I can change this
>> easily though.
>
> If "Stan users who want example models" are target audience, perhaps the data should be transformed using the tranformed
> data block.

My only concern is that the transforms are scripted, but that
can be either in R or in the Stan model. I'd then like the data
set in whatever version gets read into R in its own file so that
it can be used outside of R.

I'd prefer transforms in the Stan model if the original encoding
is the natural one ("natural" is a judgment call, of course).

>> I think the target is anyone reading ARM. I agree that we should be
>> using stan instead of lm/lmer/bugs. But I'm not sure about sim() and how
>> it works (same goes for extract()) so I'm not sure how to replace sim()
>> with extract(). If you could write up an example, that'd be great; then
>> I could just follow that.
>
> Caveat: I'm a learner :-) However...
>
> sim() gets posterior simulations of sigma and beta from a 'lm', of beta from a 'glm' or 'mer' object.
> If you look at the surce code, you see that sim() uses frequentist point estimates, their covariance and R functions
> such as mvrnorm() and rchisq() to simulate a posterior distribution.

We definitely don't need that!

But then if we recode lmer, we might wind up building something
similar on the inside.

> Stan's output _is_ a (sampled) posterior distribution. If you use uninformative priors and normal likelihood functions,
> you get similar results.

Right.

> I attach an example, ARM Figure 4.1 right.

See comments above.

I'm not the one you want commenting on the R code itself!

- Bob

[Sergio Polini]

Il 18/07/2013 00:25, Bob Carpenter ha scritto:
> Thanks for the attachment.
>
> Please keep the models separate.

Sure!
It was just a quick and dirty example for Peter about what "replacing sim() with
extract()" could mean.

Sergio

[Peter Li]

Thanks! I'm going to try and finish up to Ch 11 today, but then I need
to go back and fix up a fair amount of the graphs, so you could start
with Ch12 if you want.

I also need to find a way to plot multiple graphs side by side/above
each other. I think it is possible to do it using just ggplot2
(https://code.google.com/p/gridextra/wiki/arrangeGrob) but I need to
look into this more.

If you have any other suggestions let me know!


Peter

On 7/17/13 6:39 PM, Sergio Polini wrote:
> Il 18/07/2013 00:25, Bob Carpenter ha scritto:

>> Thanks for the attachment.
>>
>> Please keep the models separate.
>

> Sure!
> It was just an example for Peter about what "replacing sim() with
> extract()" could mean.
>
> Sergio
>

[Sergio Polini]

Il 19/07/2013 18:00, Peter Li ha scritto:
> Thanks! I'm going to try and finish up to Ch 11 today, but then I need
> to go back and fix up a fair amount of the graphs, so you could start
> with Ch12 if you want.

Ch12? I'ld say that there is no need to worry about Ch12-15, because
they are 'bayesized' in Ch16-17, so what is needed is a translation from
BUGS to Stan of the models in Ch16-17.
Am I missing anything?

> I also need to find a way to plot multiple graphs side by side/above
> each other. I think it is possible to do it using just ggplot2
> (https://code.google.com/p/gridextra/wiki/arrangeGrob) but I need to
> look into this more.

Please, could you make an example? Which figures are you thinking of?
Perhaps 13.1 or 13.3?

Thanks

Sergio

[Andrew Gelman]

The Bugs models in ARM need to be rewritten for Stan. We need to put in the Matt trick, also we might need to change the priors for the variance matrices in the varying-intercept varying-slope models. And for the logistic models we should use bernoulli_logit or whatever it's called.
A

[Bob Carpenter]

>> Il 19/07/2013 18:00, Peter Li ha scritto:

>>> I also need to find a way to plot multiple graphs side by side/above
>>> each other. I think it is possible to do it using just ggplot2
>>> (https://code.google.com/p/gridextra/wiki/arrangeGrob) but I need to
>>> look into this more.

If the multiple plots are related, you can use the facets
in ggplot2. The trick's in converting everything to the right
kind of data frame.

- Bob

[Peter Li]

Ah, in that case, I need to look more closely at the chapters. As for
what figures, I was just following whatever we have data for and is
coded online (http://www.stat.columbia.edu/~gelman/arm/software/).


Peter

[Sergio Polini]

Il 19/07/2013 18:00, Peter Li ha scritto:

> I also need to find a way to plot multiple graphs side by side/above
> each other. I think it is possible to do it using just ggplot2
> (https://code.google.com/p/gridextra/wiki/arrangeGrob) but I need to
> look into this more.

Il 20/07/2013 22:21, Bob Carpenter ha scritto:>
> If the multiple plots are related, you can use the facets
> in ggplot2. The trick's in converting everything to the right
> kind of data frame.

I think that Bob is wright.
I attach an example: Figure 13.1.

Sergio

[Tim Triche, Jr.]

the 'reshape2' package is terribly handy for that sort of thing (re-using facets instead of multiple grobs).

'melt', 'acast', and 'dcast' in particular are sort of skeleton keys for such manipulations. 

the 'GGally' and 'gridExtra' packages can also be useful for these purposes.

hope this helps.  I've used 'reshape2' quite a bit in one of my packages and it's been widely applicable, for me.

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

For more options, visit https://groups.google.com/groups/opt_out.

[Peter Li]

That's exactly what we need (I think). I'm not sure if all of the plots we need to make will be like this one in the sense that they may not all be plotting the same graph several times with different data. In any case, I'll follow your example for whatever I can!

Sergio

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

[Peter Li]

The description for sim() in the ARM package online is

"This generic function gets posterior simulations of sigma and beta from a lm object, or simulations
of beta from a glm object, or simulations of beta from a mer object"

So is the idea to simulate normal_rng with mean equal to whatever the coefficients are from running NUTS and sigma equal to the average step size?



Peter

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

[Peter Li]

Actually, can't I just pull random generations from the samples themselves, since the simulations themselves are posterior simulations?

[Bob Carpenter]

On 7/22/13 9:18 AM, Peter Li wrote:
> The description for sim() in the ARM package online is
>
> "This generic function gets posterior simulations of sigma and beta from a lm object, or simulations
> of beta from a glm object, or simulations of beta from a mer object"
>
> So is the idea to simulate normal_rng with mean equal to whatever the coefficients are from running NUTS and sigma equal
> to the average step size?

The idea is to simulate new data based on the parameters
fit by the model. Andrew can probably point you to the relevant
section of ARM and/or BDA where he discusses this.

You can do this in the generated quantities block.
Or you could do it in the model block.

- Bob

[Bob Carpenter]

On 7/22/13 9:34 AM, Peter Li wrote:
> Actually, can't I just pull random generations from the samples themselves, since the simulations themselves are
> posterior simulations?

Those are draws from the posterior, not simulated data.

Suppose you have data y and parameters mu and sigma for a simple
normal.

data {
...
vector[N] y;
model {
...
y ~ normal(mu,sigma);
}

then you want to generate new data y-twiddle, by

y_twiddle ~ normal(mu,sigma);

You can do this in the generated quantities block:

generated quantities {
vector[N_twiddle] y_twiddle;

for (i in 1:size(y_twiddle))
y_twiddle[i] <- normal_rng(mu,sigma);

where mu and sigma are the parameters fit by the model.

- Bob

> stan-users+unsubscribe@__googlegroups.com <mailto:stan-users%2Bunsu...@googlegroups.com>.
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

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

[Peter Li]

So in the case that I want to simulate the coefficient of a linear model
y ~ x, I would have something like

generated quantities {
real b1_twiddle;
real b2_twiddle;
b1_twiddle <- normal_rng(beta[1],sigma);
b2_twiddle <- normal_rng(beta[2],sigma);
}

But when I run this and try to plot each simulation as lines, they're
not very close to the data (When I run sim() from ARM, the lines go
right through the data).

[Bob Carpenter]

On 7/23/13 9:08 AM, Peter Li wrote:
> So in the case that I want to simulate the coefficient of a linear model y ~ x, I would have something like



> generated quantities {
> real b1_twiddle;
> real b2_twiddle;
> b1_twiddle <- normal_rng(beta[1],sigma);
> b2_twiddle <- normal_rng(beta[2],sigma);
> }
>
> But when I run this and try to plot each simulation as lines, they're not very close to the data (When I run sim() from
> ARM, the lines go right through the data).

No, this is simulating new data.

If all you want are samples from the posterior, you can just
use the samples returned by Stan.

- Bob

[Avi Feller]

Hi all,

Apologies for jumping onto an old chain. I was implementing a model like those in ARM Ch. 17.1 (especially those that model the correlation between intercepts and slopes) when I saw this note that some of this code might already be off-the-shelf. 

I'm having a hard time following the rest of the chain---does anyone have Stan code that they can share for ARM ch. 17 with all the latest tips/tricks/etc.? From what I can gather, the examples in stan/src/models/ARM are currently through ch. 11 (and I'm still grateful for that!).

Many thanks,

Avi

On Tuesday, July 16, 2013 2:15:03 PM UTC-4, Peter Li wrote:

The code you've attached looks great. I'm actually working on writing
the ARM models in Stan too. What I've done so far is in
https://github.com/stan-dev/stan/tree/feature/ARM. I need to do a bit of
reorganizing, but I'll mostly be following a similar format/organization
as you.

Alexij is also working on the ARM models, and from what I understand,
he's already done the majority of ~ch16-20 and I've got a fair amount of
the first 12 chapters done (without the R code). In any case, if you
want to work together on this let me know and we can coordinate who's
doing what so we don't end up writing up the same models several times. :-)


Peter

On 7/14/13 12:30 PM, sp_r...@yahoo.it wrote:
> I've read a lot of messages and noticed that:
> a) you use stopifnot(require(rstan)) instead of library(rstan);
> b) another developer (I can't remember his name...) appends .RData to
> saved compiled models.
> Moreover,
> c) as far I can understand, when using Stan from the command line the
> data file can contain more variables than the ones declared in the
> data block, so I enclose just one data file, even if only a subset of
> its variables are meant to be read by a model;
> d) remembering your words (my R code should run out of the box) I've
> rewritten my scripts; I'm used to execute a script line by line (or
> region by region) in Emacs+ESS, but you can now type source(<script
> name>) and look at the numerical and graphical results.
>
> So I've removed my previous attachment (GH2007_Ch.3.tar.gz) and attach
> a new one: GH2007_Ch.3_0.2.tar.gz
>
> Best wishes,
>
> Sergio

[Bob Carpenter]

I think Alexej and Peter may have some examples.
If not, they'll get there.

Is there a specific question we could answer now?

We should add a covariance prior example to the manual
regression chapter (in the hierarchical regression section).
Peter --- could you do that when you get to the relevant models
in ARM? We want to illustrate scaling an lkj_correlation()
prior.

- Bob

On 7/25/13 12:31 PM, Avi Feller wrote:
> Hi all,
>
> Apologies for jumping onto an old chain. I was implementing a model like those in ARM Ch. 17.1 (especially those that
> model the correlation between intercepts and slopes) when I saw this note that some of this code might already be
> off-the-shelf.
>
> I'm having a hard time following the rest of the chain---does anyone have Stan code that they can share for ARM ch. 17
> with all the latest tips/tricks/etc.? From what I can gather, the examples in stan/src/models/ARM are currently through
> ch. 11 (and I'm still grateful for that!).
>
> Many thanks,
> Avi
>
>
> On Tuesday, July 16, 2013 2:15:03 PM UTC-4, Peter Li wrote:
>
> The code you've attached looks great. I'm actually working on writing
> the ARM models in Stan too. What I've done so far is in

> https://github.com/stan-dev/stan/tree/feature/ARM <https://github.com/stan-dev/stan/tree/feature/ARM>. I need to do

> a bit of
> reorganizing, but I'll mostly be following a similar format/organization
> as you.
>
> Alexij is also working on the ARM models, and from what I understand,
> he's already done the majority of ~ch16-20 and I've got a fair amount of
> the first 12 chapters done (without the R code). In any case, if you
> want to work together on this let me know and we can coordinate who's
> doing what so we don't end up writing up the same models several times. :-)
>
>
> Peter
>

> On 7/14/13 12:30 PM, sp_r...@yahoo.it <javascript:> wrote:
> > I've read a lot of messages and noticed that:
> > a) you use stopifnot(require(rstan)) instead of library(rstan);
> > b) another developer (I can't remember his name...) appends .RData to
> > saved compiled models.
> > Moreover,
> > c) as far I can understand, when using Stan from the command line the
> > data file can contain more variables than the ones declared in the
> > data block, so I enclose just one data file, even if only a subset of
> > its variables are meant to be read by a model;
> > d) remembering your words (my R code should run out of the box) I've
> > rewritten my scripts; I'm used to execute a script line by line (or
> > region by region) in Emacs+ESS, but you can now type source(<script
> > name>) and look at the numerical and graphical results.
> >
> > So I've removed my previous attachment (GH2007_Ch.3.tar.gz) and attach
> > a new one: GH2007_Ch.3_0.2.tar.gz
> >
> > Best wishes,
> >
> > Sergio
> >
> > --
> > 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 <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Avi Feller]

Thanks for the response. Mostly, I'm looking for an example model of a hierarchical linear model that also models the correlation between varying intercepts and slopes (or, more generally, a good example with a scaled inverse-Wishart model in Stan; the example from the manual doesn't seem to get at all of this). Reading through the users group emails, I see references to all sorts of tips and tricks (I hear Matt has a trick of some sort). 

Rather than try to parse it all out on my own, I was hoping someone had a good off-the-shelf example that I could use to try to understand it all. But this is definitely not urgent---just thought I'd ask.

Thanks, as always, for the help.

Best,

Avi

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

For more options, visit https://groups.google.com/groups/opt_out.

--
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/eSYa0vTelw4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Andrew Gelman]

Hi, Avi.

You could start with a vayring-intercept varying slope model with a 2x2 matrix and no prior at all, which corresponds to a flat prior on the Cholesky parameterization, I believe (I can't remember quite what parameterization we use).  Or you could try the LKJ model which is implemented and should be explained in the stan manual.  I hope this helps!

See you

A

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

[Peter Li]

If I recall correctly, Alexej has already coded the models from Ch17. In any case, we're still working on getting everything on github and indexing all the models that we have done.


Peter

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/groups/opt_out.

--
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/eSYa0vTelw4/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/groups/opt_out.

--

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.

[Scott Baldwin]

I've created a random intercept and slope model based on "Model B"
from the longitudinal analysis textbook examples on UCLA's applied
statistics website:
http://www.ats.ucla.edu/stat/stata/examples/alda/chapter4/aldastatach4.htm

The files can be downloaded from github:
https://github.com/psychstatistics/stan_examples

I've got one other simple example on there. I am building a set of
examples for my students and for the Stan community. Just getting
started so I just have a few things going.

Anyhow, I know there are better ways to parameterize the model.
However, I have this set-up to mimic the standard frequentist
mixed-model to help the transition for my students. If anyone is
interested, I'd welcome alterations to the code that would make it
better or additional examples. Also if I have anything wrong, I'd
appreciate corrections.

Hope this is useful.

Best,
Scott

[sp_r...@yahoo.it]

I've rewritten the models, data files, and R scripts for Ch3-4 following the guidelines that I've guessed by reading your advices:
a) models in separate files;
b) data preperation (removing cases with missing values etc.) before stan_rdump()ing;
c) data transformation in the transformed data block;
d) each variable declared in a model is defined in the data file and has the same name;
e) standardization to get faster sampling (almost always, not needed when the sample size is very small);
etc.

The guidelines are detailed in a README file.
If they are agreed upon, I'll go on. Otherwise I'll adapt myself ;-)

Thanks

Sergio

[Peter Li]

These look really great. I'll need to adapt what I've done (and clean up
my code a bit) to mirror your guidelines. Thanks for putting this together!


Peter

P.S. Alexij said he'll try to have Ch 16-18 on github by the end of the
week.

[Sergio Polini]

I attach a few suggestions about Ch5. The main changes are:

1. �5.1-5.7

I've centered the variables even when not required by the model in order
to improve both speed and effective sample size. As far as I can
understand, this happens because centering reduces posterior correlation
of the parameters. Parameter values are recovered in the generated
quantities block.
However:
-- I don't know if standardization would be better (I did not try);
-- I can't understand why the effective sample size gets larger for
parameters but smaller for lp__.

2. �5.1,5.2

a) I've splitted the nes dataset.

b) I've used geom_line() instead of geom_smooth(method = "glm") because
geom_smooth() do.call()s the glm() function, so using it when plotting
Stan results is somewhat like cheating ;-)

4. �5.4-5.6

I've annotated the figures.

5. �5.6

I've generated the predicted (fitted) values in the generated quantities
block.
But I'm not sure that it's a good idea.

HTH

Sergio

[Sergio Polini]

A few suggestions about Ch9.

Sergio

[Peter Li]

Great! I'll look at both changes for Ch5 and Ch9 and update whatever
files are necessary. Thanks again for doing this. And for anyone else
that is interested in the Stan code for ARM, you can find an index of
what has been completed and on github at
https://github.com/stan-dev/stan/wiki/ARM-Models.



Peter

[Bob Carpenter]

Wow! This looks awesome. I hope the indexing
wasn't too onerous --- I really like having it both
ways.

I do have one minor suggestion: rather than having the
link anchor text for all the Kid IQ model variants
be "kid_iq", the anchor text should match the target file
name, like "kid_iq_c" and "kid_iq_interaction".

- Bob

[Bob Carpenter]

On 8/9/13 3:54 AM, Sergio Polini wrote:
> I attach a few suggestions about Ch5. The main changes are:
>
> 1. �5.1-5.7
>
> I've centered the variables even when not required by the model in order to improve both speed and effective sample
> size.

Do we want to do this or have the model examples follow
the text?

> As far as I can understand, this happens because centering reduces posterior correlation of the parameters.

Linear transformation like centering and scaling won't
affect correlation. But they will affect the intercept,
which has to absorb the slack from offsets. Look
at what happens to predictors under centering:

alpha * x[i] + beta

== alpha * (x[i] - mean(x)) + (beta + alpha * mean(x))

The coefficient is the same, but the new intercept is offset
by alpha * mean(x).

> Parameter values are recovered in the generated quantities block.
> However:
> -- I don't know if standardization would be better (I did not try);

If the variables have different scales, then yes. We do adaptation,
but it's approximate and it's faster to just start with everything on the
same unit scale. This does affect the meaning of the priors, though.

> -- I can't understand why the effective sample size gets larger for parameters but smaller for lp__.

We don't really care about ESS for lp__ so much. It's a
non-linear function of the parameters (the log probability
function to be exact).

> 2. �5.1,5.2
>
> a) I've splitted the nes dataset.
>
> b) I've used geom_line() instead of geom_smooth(method = "glm") because geom_smooth() do.call()s the glm() function, so
> using it when plotting Stan results is somewhat like cheating ;-)

What are you plotting? An alternative would be to use
something like loess.

> 4. �5.4-5.6
>
> I've annotated the figures.
>
> 5. �5.6
>
> I've generated the predicted (fitted) values in the generated quantities block.
> But I'm not sure that it's a good idea.

Yes, that's the best place to do it both for clarity of the model
and for efficiency.

Thanks for all the work on this. As I wrote in the previous message,
it all looks great.

- Bob

[Andrew Gelman]

My suggestion is:

Matt trick: Yes.

Centering predictors (when they were not centered in the original model): No. Unless Stan is running really slow with non-centered. But if it runs fast with non-centered I'd prefer to keep the predictors as is.

A

[Sergio Polini]

Il 09/08/2013 19:28, Bob Carpenter ha scritto:
> On 8/9/13 3:54 AM, Sergio Polini wrote:
>> b) I've used geom_line() instead of geom_smooth(method = "glm")
>> because geom_smooth() do.call()s the glm() function, so
>> using it when plotting Stan results is somewhat like cheating ;-)
>
> What are you plotting? An alternative would be to use
> something like loess.

I was plotting figures like 5.1 (a).
However you always need a set of (x,y) points.
If you use geom_smooth(method="glm"), the y values are calculated by
predict.glm(model, ...) where model is an object returned by glm() -- at
least, this is what I guess by looking at ggplot2/R/stat-smooth.r.
The y values predicted by Stan are very similar if you use vague priors,
but could be different.
This is why I've used:

beta.mean <- colMeans(extract(<stanfit>, "beta")$beta)
y <- 1 / (1 + exp(- beta.mean[1] - beta.mean[2] * x))

and then geom_line(), something like the example in the geom_smooth()
help page.

> Thanks for all the work on this. As I wrote in the previous message,
> it all looks great.

I thank you, because I'm learning a lot!

Sergio

[sp_r...@yahoo.it]

Ok. I've rewritten and attach the Stan code for Ch3-5,9 leaving centering and standardization only when required by the models.

Sergio

[Peter Li]

Just finished making the fixes. I'll update the other pages once all the
models are in (Ch16-18 still missing).


Peter

[Peter Li]

Great! I just pushed your fixes to the github repo.


Peter

[sophie]

Hi Peter,

Did you also upload the R codes and codes for ggplot2 for Gelman and Hill 03 book? I only see the Stan codes on github. Thanks!

Best,

Yajuan

[Peter Li]

Yes. They are all on github for all chapters except 16-18 [those should be on soon]. If you go to https://github.com/stan-dev/stan/wiki/ARM-Models-Sorted-by-Chapter, the hyperlinks for each section (i.e. Ch 2.3 CI) will lead you to the part of the github repo with the R code/ggplot code for that section. Some sections have been omitted due to the lack of data.


Peter

[sp_r...@yahoo.it]

I attach what I've done about Ch.16.

ARM> 16.3_Fitting & understanding a varying-intercept mlm using R and Bugs.R
+ 16.4_Step by step through a Bugs model.R
+ 16.8_Principles of modeling in Bugs.R

I've replaced "Bugs" with "Stan", "R and Bugs" with "RStan".
I hope that the content (some topics discussed below) still makes sense.

ARM> 16.4_Step by step through a Bugs model.R

I've not converted:
radon.1.noburnin <- bugs (..., n.burnin=0)
because if I try:
sf <- sampling(..., warmup = 0) # or warmup = 1
I get no convergence.
I've just used extract(..., permuted = FALSE, inc_warmup = TRUE).

ARM> 16.4_Bugs_codes.bug

The Stan code is in radon.1.stan

ARM> 16.5_Adding individual & group-level predictors.R
+ 16.5_Bugs_codes.bug
+ 16.5_vector.matrix.notation.intercepts.bug

The first file looks like "16.4_Step by step through a Bugs model.R" without the code for figure 16.2 and with a few more lines of R code.
I've named "16.5_AddingIndividualAndGroup-levelPredictors.R" a file containing R code to fit:
- radon.pooling.stan (1st model in "16.5_Bugs_codes.bug")
- radon.nopooling.stan (2nd model in "16.5_Bugs_codes.bug")
- radon.2.stan (last model in "16.5_Bugs_codes.bug")
The data for the other Bugs models are missing.

ARM> 16.6_Predictions.R

The data are augmented with records containing missing values, but Stan doesn't support arrays of known and missing data.
So I've not rewritten radon.2a.bug to a unique radon.2a.stan, but I've written:
a) radon.2a.stan: a new unit in an existing group;
b) radon.2b.stan: a new unit in a new group.
In both cases, y_tilde is generated in the generated quatities block.

ARM> 16.8_Principles\ of\ modeling\ in\ Bugs.R

radon.3.bug is just alike radon.2.bug but for the second commented line.  Moreover, if one wants to specify the values of parameters when using Stan, the unmodeled parameters must be declared in the data or transformed data block.
So I've used radon.2.stan (rewriting of radon.2.bug) for the first sampling, then a new radon.3.stan with a transformed data block:

transformed data {
  real<lower=0> sigma_y;
  real<lower=0> sigma_a;
  sigma_y <- .7;
  sigma_a <- .4;
}

ARM> 16.9_Practical issues of implementation.R + 16.9_Bugs_codes.bug
+ 16.10_Open-ended modeling in Bugs.R

16.9_* contain some code for a schools example, which is in §16.10.
16_10 contains examples of inits() functions, which are in §16.9.
I think that titles and contents should be switched.
As to attached file:
a) I do not include a 16.9 file, because the inits() functions refer to hypothetical data;
b) I do not include a 16.10 file, because as far as I can understand the schools examples is again a hypothetical one, neither the Bugs ranking schools example, nor the eigth-schools one.

HTH

Sergio

[Bob Carpenter]

Thanks --- is this in version control somewhere
I can check out?

More inline.

On 8/25/13 5:42 AM, sp_r...@yahoo.it wrote:
> I attach what I've done about Ch.16.
>
> ARM> 16.3_Fitting & understanding a varying-intercept mlm using R and Bugs.R
> + 16.4_Step by step through a Bugs model.R
> + 16.8_Principles of modeling in Bugs.R
>
> I've replaced "Bugs" with "Stan", "R and Bugs" with "RStan".
> I hope that the content (some topics discussed below) still makes sense.
>
> ARM> 16.4_Step by step through a Bugs model.R
>
> I've not converted:
> radon.1.noburnin <- bugs (..., n.burnin=0)
> because if I try:
> sf <- sampling(..., warmup = 0) # or warmup = 1
> I get no convergence.

You mean in Stan? We don't want to run without warmup.
And Andrew wants the terminology changed form "burnin"
to "warmup" --- it's a better physical analogy.

> I've just used extract(..., permuted = FALSE, inc_warmup = TRUE).
>
> ARM> 16.4_Bugs_codes.bug
>
> The Stan code is in radon.1.stan
>
> ARM> 16.5_Adding individual & group-level predictors.R
> + 16.5_Bugs_codes.bug
> + 16.5_vector.matrix.notation.intercepts.bug
>
> The first file looks like "16.4_Step by step through a Bugs model.R" without the code for figure 16.2 and with a few
> more lines of R code.
> I've named "16.5_AddingIndividualAndGroup-levelPredictors.R" a file containing R code to fit:
> - radon.pooling.stan (1st model in "16.5_Bugs_codes.bug")
> - radon.nopooling.stan (2nd model in "16.5_Bugs_codes.bug")
> - radon.2.stan (last model in "16.5_Bugs_codes.bug")

I'd rather have individual R files for each Stan model
that read in data and fit and produce outputs.

You could then collect these into bigger chunks, but
I don't want to have to understand a bigger piece of script
than I need to in order to fit a model.

Other R users may have different opinions on how to do this.

> The data for the other Bugs models are missing.

Missing as in you can't find it, or missing as in missing
data models?

> ARM> 16.6_Predictions.R
>
> The data are augmented with records containing missing values, but Stan doesn't support arrays of known and missing data.

Not directly as input, but check out the manual chapter on
missing data to see how to encode missing data problems.

> So I've not rewritten radon.2a.bug to a unique radon.2a.stan, but I've written:
> a) radon.2a.stan: a new unit in an existing group;
> b) radon.2b.stan: a new unit in a new group.
> In both cases, y_tilde is generated in the generated quatities block.

Sounds good.

> ARM> 16.8_Principles\ of\ modeling\ in\ Bugs.R
>
> radon.3.bug is just alike radon.2.bug but for the second commented line. Moreover, if one wants to specify the values
> of parameters when using Stan, the unmodeled parameters must be declared in the data or transformed data block.

By unmodeled parameter, you just mean a constant?
Like if I did:

theta ~ beta(alpha,beta);

and declared alpha and beta as data and gave them
fixed values?

> So I've used radon.2.stan (rewriting of radon.2.bug) for the first sampling, then a new radon.3.stan with a transformed
> data block:
>
> transformed data {
> real<lower=0> sigma_y;
> real<lower=0> sigma_a;
> sigma_y <- .7;
> sigma_a <- .4;
> }

I like this like this so that you can see the assumptions
in the model. The other thing to do would be to just remove
sigma_y and use 0.7 in its place, but the naming can help
with understanding if it gets used in more than one place.

> ARM> 16.9_Practical issues of implementation.R + 16.9_Bugs_codes.bug
> + 16.10_Open-ended modeling in Bugs.R
>
> 16.9_* contain some code for a schools example, which is in �16.10.
> 16_10 contains examples of inits() functions, which are in �16.9.
> I think that titles and contents should be switched.
> As to attached file:
> a) I do not include a 16.9 file, because the inits() functions refer to hypothetical data;
> b) I do not include a 16.10 file, because as far as I can understand the schools examples is again a hypothetical one,
> neither the Bugs ranking schools example, nor the eigth-schools one.

Not sure what you mean by "hypothetical data". Did you
mean simulated? Eight schools is included in our install guide.
And also discussed by Andrew in the appendix to the next edition
of BDA.

- Bob

[Sergio Polini]

Il 25/08/2013 22:32, Bob Carpenter ha scritto:
> Thanks --- is this in version control somewhere
> I can check out?

Sorry, I can't understand. Were you asking Peter?

>> ARM> 16.4_Step by step through a Bugs model.R
>>
>> I've not converted:
>> radon.1.noburnin <- bugs (..., n.burnin=0)
>> because if I try:
>> sf <- sampling(..., warmup = 0) # or warmup = 1
>> I get no convergence.
>
> You mean in Stan? We don't want to run without warmup.
> And Andrew wants the terminology changed form "burnin"
> to "warmup" --- it's a better physical analogy.

Yes, in Stan. No problem, however. Using inc_warmup I was able to plot
the first 200 iterations, as required by the ARM code.

>> I've named "16.5_AddingIndividualAndGroup-levelPredictors.R" a file
>> containing R code to fit:
>> - radon.pooling.stan (1st model in "16.5_Bugs_codes.bug")
>> - radon.nopooling.stan (2nd model in "16.5_Bugs_codes.bug")
>> - radon.2.stan (last model in "16.5_Bugs_codes.bug")
>
> I'd rather have individual R files for each Stan model
> that read in data and fit and produce outputs.
>
> You could then collect these into bigger chunks, but
> I don't want to have to understand a bigger piece of script
> than I need to in order to fit a model.
>
> Other R users may have different opinions on how to do this.

The R files can be rewritten (splitted), but they often contain several
models because so do the ones in
http://www.stat.columbia.edu/~gelman/arm/examples/Book_Codes.zip

>> The data for the other Bugs models are missing.
>
> Missing as in you can't find it, or missing as in missing
> data models?

Missing because there is not a full example, but just an outline.
For instance, in 16.5_Bugs_codes.bug there is a model (Classical
regression with multiple predictors) with two input variables, x and
winter, but the second one is hypothetical, it is not in the radon data.

>> ARM> 16.8_Principles\ of\ modeling\ in\ Bugs.R
>>
>> radon.3.bug is just alike radon.2.bug but for the second commented
>> line. Moreover, if one wants to specify the values
>> of parameters when using Stan, the unmodeled parameters must be
>> declared in the data or transformed data block.
>
> By unmodeled parameter, you just mean a constant?
> Like if I did:
>
> theta ~ beta(alpha,beta);
>
> and declared alpha and beta as data and gave them
> fixed values?

Yes. I mean "const. unmodeled param" as in the manual, page 181.

>> ARM> 16.9_Practical issues of implementation.R + 16.9_Bugs_codes.bug
>> + 16.10_Open-ended modeling in Bugs.R
>>
>> 16.9_* contain some code for a schools example, which is in �16.10.
>> 16_10 contains examples of inits() functions, which are in �16.9.
>> I think that titles and contents should be switched.
>> As to attached file:
>> a) I do not include a 16.9 file, because the inits() functions refer
>> to hypothetical data;
>> b) I do not include a 16.10 file, because as far as I can understand
>> the schools examples is again a hypothetical one,
>> neither the Bugs ranking schools example, nor the eigth-schools one.
>
> Not sure what you mean by "hypothetical data". Did you
> mean simulated? Eight schools is included in our install guide.
> And also discussed by Andrew in the appendix to the next edition
> of BDA.

The schools example in �16.10 (ARM, page 370) is about "a hypothetical
study of a new teaching method applied in J different classrooms
containing a total of n students. Our data for this example will be the
treatment indicator T (defined at the school level) and, for each
student, a pre-treatment assessment, x (on a 1-10 scale, say) and a
post-treatment test score, y (on a 0-100 scale)".
Neither schools (Bugs) nor eight-schools contain such data.

Thanks for your feedback.

Sergio

[Peter Li]

Sergio, are you on github? If so, I think we can work something out so
you can directly add your code to feature/ARM. If not, I can add the
files. I can also take a look at these files this week. Also, I believe
Alexej has Ch 16 codes, but hasn't put them on github yet.

[Sergio Polini]

Il 26/08/2013 04:57, Peter Li ha scritto:
> Sergio, are you on github? If so, I think we can work something out so
> you can directly add your code to feature/ARM. If not, I can add the
> files. I can also take a look at these files this week.

No, I'm not on github.

[sp_r...@yahoo.it]

I attach data and code for Ch18, i.e. for "18.3_Bayes for classical & multilevel regression.R" and "18.5_Likelihood, Bayes & Gibbs for censored data.R."
I have not rewritten "18.4_Gibbs sampler for mlm.R" nor "18.7_Programming Gibbs & Metropolis in R.R" because:
1) I do not know how to code HMC/NUTS in R. I don't even know if it's feasible or makes sense.
2) I can't find the data for the social network example.

Sergio

[Bob Carpenter]

On 8/28/13 5:37 AM, sp_r...@yahoo.it wrote:
> I attach data and code for Ch18, i.e. for "18.3_Bayes for classical & multilevel regression.R" and "18.5_Likelihood,
> Bayes & Gibbs for censored data.R."
> I have not rewritten "18.4_Gibbs sampler for mlm.R" nor "18.7_Programming Gibbs & Metropolis in R.R" because:

Yup -- those don't need to be changed.

It's easy enough to code up HMC directly. NUTS is quite
a bit more work. Radford Neal probably has some HMC in R
code we could reference:

http://www.cs.toronto.edu/~radford/software-online.html

http://www.cs.toronto.edu/~radford/ham-mcmc.software.html

It looks like he's unifying his interfaces and rollling things into
his GRIMS package:

http://www.cs.toronto.edu/~radford/GRIMS.html

> 1) I do not know how to code HMC/NUTS in R. I don't even know if it's feasible or makes sense.
> 2) I can't find the data for the social network example.

If Andrew doesn't answer this, send him mail directly with
a relevant subject.

- Bob

[Andrew Gelman]

Hi--direct R code of HMC is in appendix C of BDA3. No need to have it in ARM, it's too complicated for that book.
A

[Sergio Polini]

I'm writing the code for Ch17, but I'm in trouble...
I'm trying to rewrite wishart1.bug (data and models attached.)
My first attempt:

data {
int<lower=0> N;
int<lower=0> J;
real y[N];
real x[N];
int county[N];
}
transformed data {
cov_matrix[2] Scale;
int df;
for (i in 1:2)
for (j in 1:2)
if (i == j) Scale[i,j] <- 1.0; else Scale[i,j] <- 0.0;
df <- 3;
}
parameters {
real<lower=0> xi_a;
real<lower=0> xi_b;
vector[2] B_raw[J];
vector[2] mu_B_raw;
real<lower=0> sigma_y;
cov_matrix[2] Sigma2_B_raw;
}
transformed parameters {
vector[J] a;
vector[J] b;
real mu_a;
real mu_b;
real sigma_a;
real sigma_b;
for (j in 1:J) {
a[j] <- xi_a * B_raw[j,1];
b[j] <- xi_b * B_raw[j,2];
}
mu_a <- xi_a * mu_B_raw[1];
mu_b <- xi_b * mu_B_raw[2];
sigma_a <- xi_a * sqrt(Sigma2_B_raw[1,1]);
sigma_b <- xi_b * sqrt(Sigma2_B_raw[2,2]);
}
model {
Sigma2_B_raw ~ inv_wishart(df, Scale);
for (j in 1:J)
B_raw[j] ~ multi_normal(mu_B_raw, Sigma2_B_raw);
for (i in 1:N)
y[i] ~ normal(a[county[i]] + b[county[i]] * x[i], sigma_y);
}
generated quantities {
real<lower=-1,upper=1> rho;
rho <- Sigma2_B_raw[1,2] / sqrt(Sigma2_B_raw[1,1] * Sigma2_B_raw[2,2]);
}

It works, but even after 20,000 iterations:
a) effective sample size is very low: about 50 for x_a and B_raw[*,1], 7
for x_b, 3 or 4 for B_raw[*,2], a few hundreds for most of the others, 8
for sigma_b;
b) Rhat is greater than 1.7 for B_raw[*,2], 1.2 for sigma_b.

I must say that I don't understand well how statements in the
transformed parameters block work.
For instance, may I define a[j] as xi_a * B_raw[j,1] when B_raw is not
yet defined, because it is modeled in the model block?
Yes, I think, but I've tried another way: no transformed parameters
block, "a" and "b" replaced with "xi_a * B_raw" in the model block, the
derived quantities in the generated quantities block:

...

model {
Sigma2_B_raw ~ inv_wishart(df, Scale);
for (j in 1:J)
B_raw[j] ~ multi_normal(mu_B_raw, Sigma2_B_raw);
for (i in 1:N)
y[i] ~ normal(xi_a * B_raw[county[i],1] + xi_b * B_raw[county[i],2]
* x[i], sigma_y);
}
generated quantities {
vector[J] a;
vector[J] b;
real mu_a;
real mu_b;
real sigma_a;
real sigma_b;
real<lower=-1,upper=1> rho;
for (j in 1:J) {
a[j] <- xi_a * B_raw[j,1];
b[j] <- xi_b * B_raw[j,2];
}
mu_a <- xi_a * mu_B_raw[1];
mu_b <- xi_b * mu_B_raw[2];
sigma_a <- xi_a * sqrt(Sigma2_B_raw[1,1]);
sigma_b <- xi_b * sqrt(Sigma2_B_raw[2,2]);
rho <- Sigma2_B_raw[1,2] / sqrt(Sigma2_B_raw[1,1] * Sigma2_B_raw[2,2]);
}

Result (after 20,000 iterations): lower ESS, higher Rhat!

What about reparametrization?
It should be simple, because the Scale matrix is an identity matrix, S =
I2, so:
a) Cholesky factor: L = I2;
b) A = [ sqrt(c1), 0 ; z, sqrt(c2) ];
c) inv(W) = inv(L') * inv(A') * inv(A) * inv(L) = inv(A') * inv(A)
= t(inv(A)) * inv(A) ~ inverse-Wishart.
I've modified the first model (derived quantities in the transformed
parameters block):

model {
matrix[2,2] A;
matrix[2,2] A_inv_t;
z ~ normal(0, 1);
for (i in 1:2) {
c[i] ~ chi_square(df - i + 1);
A[i,i] <- sqrt(c[i]);
}
A[2,1] <- z;
A[1,2] <- 0.0;
A_inv_t <- mdivide_left_tri_low(A, I2); // I2 = identity matrix
A_inv_t <- A_inv_t';
for (j in 1:J)
B_raw[j] ~ multi_normal_cholesky(mu_B_raw, A_inv_t);
for (i in 1:N)
y[i] ~ normal(a[county[i]] + b[county[i]] * x[i], sigma_y);
}

Result (after 20,000 iterations):
a) c: ESS = 25 (c[1]) and 6 (c[2]), Rhat > 1.2;
b) B_raw[*,1]: ESS <= 35, Rhat > 1.15;
c) B_raw[*,2]: ESS <= 6, Rhat =~ 1.6;
d) xi_a and xi_b: ESS = 11, Rhat >= 1.24.

Moreover, I can't understand how to recover the covariance matrix, the
old Sigma2_B_raw.

I hope someone can help me ;-)

Thanks

Sergio

[Bob Carpenter]

On 9/3/13 10:22 AM, Sergio Polini wrote:
> I'm writing the code for Ch17, but I'm in trouble...

Uh oh. Let's see if we can help.

> I'm trying to rewrite wishart1.bug (data and models attached.)
> My first attempt:
>
> data {
> int<lower=0> N;
> int<lower=0> J;
> real y[N];
> real x[N];
> int county[N];
> }
> transformed data {
> cov_matrix[2] Scale;
> int df;
> for (i in 1:2)
> for (j in 1:2)
> if (i == j) Scale[i,j] <- 1.0; else Scale[i,j] <- 0.0;
> df <- 3;
> }

You can replace the loops with this:

Scale <- diag_matrix(rep_vector(1.0,2));

I like the idea of naming variables like df in the transformed
data block.

This will be faster given that it's in a loop if you first
pull out the Cholesky factor of the scale and then use
multi_normal_cholesky:

Matrix[2,2] L;
L <- cholesky_decompose(Sigma2_B_raw);
...
multi_normal_cholesky(mu_B_raw, L);

> for (i in 1:N)
> y[i] ~ normal(a[county[i]] + b[county[i]] * x[i], sigma_y);

This will be faster if you put the predictors into a vector
and then call the vectorized version of normal. Loops are fast
in Stan, but calls to probability distributions are not.

> generated quantities {
> real<lower=-1,upper=1> rho;
> rho <- Sigma2_B_raw[1,2] / sqrt(Sigma2_B_raw[1,1] * Sigma2_B_raw[2,2]);
> }

> It works, but even after 20,000 iterations:
> a) effective sample size is very low: about 50 for x_a and B_raw[*,1], 7 for x_b, 3 or 4 for B_raw[*,2], a few hundreds
> for most of the others, 8 for sigma_b;

It's a tough model. This is the kind of model where we might
want to standardize (a multi-dimensioanl Matt trick).

> b) Rhat is greater than 1.7 for B_raw[*,2], 1.2 for sigma_b.

How are you initializing? We want to run longer than this.

What is rho on the output? If it's very high (close to 1), it's going to be
problematic for sampling.

> I must say that I don't understand well how statements in the transformed parameters block work.

Have you read the manual language sections on transformed parameters?

In a nutshell, all of the data, transformed data, and parameters are
available when you compute transformed parameters, which will then be
available in the model.

> For instance, may I define a[j] as xi_a * B_raw[j,1] when B_raw is not yet defined, because it is modeled in the model
> block?

B_raw is declared in the parameters block, so it will be
defined at the point the transformed parameter block is
executed. The sampling statmeents do NOT define the variables
on their left side, they just increment the total log probability
based on variables that already have to be defined.

> Yes, I think, but I've tried another way: no transformed parameters block, "a" and "b" replaced with "xi_a * B_raw" in
> the model block, the derived quantities in the generated quantities block:

That'll work, too. If you don't want the output, move
the transformed parameters into local variables in the model,
or just get rid of them altogether.

> ...
>
> model {
> Sigma2_B_raw ~ inv_wishart(df, Scale);
> for (j in 1:J)
> B_raw[j] ~ multi_normal(mu_B_raw, Sigma2_B_raw);
> for (i in 1:N)
> y[i] ~ normal(xi_a * B_raw[county[i],1] + xi_b * B_raw[county[i],2] * x[i], sigma_y);
> }
> generated quantities {
> vector[J] a;
> vector[J] b;
> real mu_a;
> real mu_b;
> real sigma_a;
> real sigma_b;
> real<lower=-1,upper=1> rho;
> for (j in 1:J) {
> a[j] <- xi_a * B_raw[j,1];
> b[j] <- xi_b * B_raw[j,2];
> }
> mu_a <- xi_a * mu_B_raw[1];
> mu_b <- xi_b * mu_B_raw[2];
> sigma_a <- xi_a * sqrt(Sigma2_B_raw[1,1]);
> sigma_b <- xi_b * sqrt(Sigma2_B_raw[2,2]);
> rho <- Sigma2_B_raw[1,2] / sqrt(Sigma2_B_raw[1,1] * Sigma2_B_raw[2,2]);
> }
>
> Result (after 20,000 iterations): lower ESS, higher Rhat!

Probably just be randomization if it's the same model.

I'd try init=0, or better yet, try sensible initial values if
you can compute them before fitting the model.

You don't need to go this far here. Just the computation
of the Cholesky factor up front as I suggest above will save a factorization
for each n in 1:N in the calls to multi_normal, as well as a lot of
error checking.

> Moreover, I can't understand how to recover the covariance matrix, the old Sigma2_B_raw.

I didn't follow the details here --- usually you just
multiply L * L' if L is the Cholesky factor for the matrix.
But I'm guessing that's not the answer to your question.

- Bob

[Andrew Gelman]

Another option is to just ditch the scaled-inverse-Wishart and use LKJ or something else. My Sigma2_B_raw etc. notation is ugly, that's for sure, and I'd be happy to ditch it.
A

[Sergio Polini]

Il 03/09/2013 21:01, Bob Carpenter ha scritto:
>
> On 9/3/13 10:22 AM, Sergio Polini wrote:
>> I'm writing the code for Ch17, but I'm in trouble...
>
> Uh oh. Let's see if we can help.

Hi Tutor!
Seriously, I thank you a lot.

But I wish to say that Andrew Gelman ("Another option is to just ditch
the scaled-inverse-Wishart and use LKJ or something else") is likely
right ;-)
I've looked around and found this discussion:

https://groups.google.com/forum/#!msg/stan-dev/21BZjOp1IAg/tOmUhMLKzmwJ

and a message by Ben Goodrich: "The {scaled, inverse} Wishart priors are
overly popular due to their computational convenience in Gibbs samplers,
but I suspect they are computationally inconvenient for HMC due to the
induced high dependence among parameters. Maybe NUTS 3 will make this
feasible, but it should already be feasible in NUTS 1 with an LKJ prior
and eta >= 1."

So replacing BUGS+Scaled-inverse-Wishart with Stan+LKJ could make sense
even when rewriting the ARM code...

I've tried to keep things as simple as possibile:
> library(MASS)
> N <- 1000 # or 10000 or...
> Sigma <- matrix(c(10,3,3,2), 2, 2)
> y <- mvrnorm(n=N, rep(0, 2), Sigma)

but all my Stan+Scaled-inverse-Wishart models were disappointing: slow
sampling, low ESS and high Rhat.
Even when I applied your suggestions, even the multidimensional Matth trick.

So I've tried:

data {
int N;
int K;
vector[K] y[N];
}
parameters {
vector[K] mu;
vector<lower=0>[K] sd;
corr_matrix[K] L;
}
transformed parameters {
matrix[K,K] S;
cov_matrix[K] Sigma;
S <- diag_matrix(sd);
Sigma <- S * L * S;
}
model {
L ~ lkj_corr(1.0);
for (n in 1:N)
y[n] ~ multi_normal(mu, Sigma);
}

No trick but LKJ.
Have I misundestood LKJ? Very likely!
But here are the results (N = 10000, iter=2000, not iter = \infty):

mean se_mean sd 50% n_eff Rhat
mu[1] 0.05 0.00 0.03 0.05 2349 1
mu[2] 0.01 0.00 0.01 0.01 2342 1
sd[1] 3.19 0.00 0.02 3.19 2256 1
sd[2] 1.42 0.00 0.01 1.42 2424 1
L[1,1] 1.00 0.00 0.00 1.00 4000 NaN
L[1,2] 0.67 0.00 0.01 0.67 2297 1
L[2,1] 0.67 0.00 0.01 0.67 2297 1
L[2,2] 1.00 0.00 0.00 1.00 1569 1
S[1,1] 3.19 0.00 0.02 3.19 2256 1
S[1,2] 0.00 0.00 0.00 0.00 4000 NaN
S[2,1] 0.00 0.00 0.00 0.00 4000 NaN
S[2,2] 1.42 0.00 0.01 1.42 2424 1
Sigma[1,1] 10.15 0.00 0.15 10.15 2257 1
Sigma[1,2] 3.06 0.00 0.06 3.06 2051 1
Sigma[2,1] 3.06 0.00 0.06 3.06 2051 1
Sigma[2,2] 2.03 0.00 0.03 2.03 2425 1
lp__ -22083.59 0.05 1.58 -22083.24 946 1

Faster sampling, reasonable means/medians, ESSs, and Rhats (as to those
NaNs, they come from a division by a zero variance when the "parameter"
is a constant, don't they?)

I'ld like such a solution.

Sergio


PS:

>> I must say that I don't understand well how statements in the
>> transformed parameters block work.
>
> Have you read the manual language sections on transformed parameters?
>
> In a nutshell, all of the data, transformed data, and parameters are
> available when you compute transformed parameters, which will then be
> available in the model.
>
>> For instance, may I define a[j] as xi_a * B_raw[j,1] when B_raw is not
>> yet defined, because it is modeled in the model
>> block?
>
> B_raw is declared in the parameters block, so it will be
> defined at the point the transformed parameter block is
> executed. The sampling statmeents do NOT define the variables
> on their left side, they just increment the total log probability
> based on variables that already have to be defined.

This is what I was missing. Thanks!

[sp_r...@yahoo.it]

I attach a file with Stan code for Ch.17.1-2.
I've replaced the scaled inverse-Wishart with implicit uniform distribution on covariance matrix (unifcov*.stan) or on correlation matrix (unifcorr*.stan).

Sergio

[Bob Carpenter]

Thanks. I'm creating a new topic to ask a followup question:

Do we want proper priors in all these models?

A uniform distribution on a correlation matrix is
proper, but not on a covariance matrix.

- Bob

[Duncan]

Well, since you've so kindly offered, Bob, I'd very much appreciate seeing the model under "17.6 Multilevel ordered categorical regression" (p. 383 of Gelman and Hill) translated into Stan!
Duncan
 

On Friday, April 5, 2013 5:09:53 AM UTC+10:30, Thomas Zumbrunn wrote:

Dear (R)Stan users

I was wondering if anyone has translated or knows of someone who has
translated all the BUGS code in

Gelman A & Hill J (2006) Data Analysis Using Regression and
Multilevel/Hierarchical Models. Cambridge University Press

to Stan code. I'm reading the book with a couple of colleagues, and we'd like
to use Stan instead of (Open|Win)BUGS. Any hints are appreciated.

Thanks,
/thomas