feedback on hierarchical model tutorial

161 views
Skip to first unread message

Bob Carpenter

unread,
Jan 15, 2016, 4:26:21 PM1/15/16
to stan...@googlegroups.com
Before releasing to our users, I'd like to get some feedback
on a knitr doc I wrote as an intro to hierarchical modeling
for repeated binary trials. It lives here:

stan-dev/example-models/knitr/efron-morris/pool-no-pool.Rmd

As usual, I had a hard time judging what to explain and what
to take for granted. I tried to focus on Stan. Any tips on
making the RStan better would be appreciated. Andrew, now's the
time to jump in on ggplot2 defaults! Or I can ask you in person
and we can create a theme you can live with.

I'm happy to send the 2.8MB HTML file it produced (encapsulated)
if you'd rather not run it yourself. (If you do want to run it,
you only need:

> rmarkdown::render("pool-no-pool.Rmd")

to get the basic HTML. I couldn't get it to keep the figures
separate, even trying what Ben suggested (I'm probably doing it
wrong).

Here's the abstract:

This note illustrates the effects on posterior inference of pooling
data (aka sharing strength) across items for repeated binary trial
data. It provides Stan models and R code to fit and check predictive
models for three situations: (a) complete pooling, which assumes each
item is the same, (b) no pooling, which assumes the items are
unrelated, and (c) partial pooling, where the similarity among the
items is estimated. We consider two hierarchical models to estimate
the partial pooling, one with a beta prior on chance of success and
another with a normal prior on the log odds of success. The note
explains with working examples how to (i) fit models in RStan and plot
the results in R ussing ggplot2, (ii) estimate event probabilities,
(iii) evaluate posterior predictive densities to compare model
predictions, (iv) rank items by chance of success, (v) perform
multiple comparisons in several settings, (vi) replicate new data for
posterior p-values, and (vii) perform graphical posterior predictive
checks.

- Bob

Aki Vehtari

unread,
Jan 16, 2016, 9:53:42 AM1/16/16
to stan development mailing list
On Friday, January 15, 2016 at 11:26:21 PM UTC+2, Bob Carpenter wrote:
> I'm happy to send the 2.8MB HTML file it produced (encapsulated)
> if you'd rather not run it yourself. (If you do want to run it,
> you only need:
>
> > rmarkdown::render("pool-no-pool.Rmd")

I got an error
quitting from lines 696-706 (pool-no-pool.Rmd)
Error in eval(expr, envir, enclos) : object 'ss_hier_logit' not found

I'm new to R (started learning a week ago), so maybe it's my fault.

You can also send me the HTML file.

Aki

Bob Carpenter

unread,
Jan 16, 2016, 1:07:21 PM1/16/16
to stan...@googlegroups.com
Thanks. That was my fault. R keeps things in the environment,
so I added a use of ss_hier_logit above where I defined it and
it worked just fine locally. I should've tried it from scratch
before sending it out.

I'll fix the issue, and I'll send you the HTML. The PDF's
just too huge because it uses PDF figures (probably configurable
if I could somehow penetrate the knitr doc).

- Bob
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Ben Goodrich

unread,
Jan 16, 2016, 1:55:38 PM1/16/16
to stan development mailing list
On Saturday, January 16, 2016 at 1:07:21 PM UTC-5, Bob Carpenter wrote:
I'll fix the issue, and I'll send you the HTML.  The PDF's
just too huge because it uses PDF figures (probably configurable
if I could somehow penetrate the knitr doc).

Yihui helped us reduce the size of the rstanarm HTML vignettes

https://github.com/stan-dev/rstanarm/pull/53

Ben





Bob Carpenter

unread,
Jan 16, 2016, 5:30:21 PM1/16/16
to stan...@googlegroups.com
I'm not sure what I'm looking at. Which piece of it
reduces size relative to the defaults?

I saw this file

vignettes/children/SETTINGS-knitr.txt

with this in it:

```{r, SETTINGS-knitr, include=FALSE}
-stopifnot(require(svglite))
stopifnot(require(knitr))
opts_chunk$set(
comment=NA, message = FALSE, warning = FALSE,
- dev = 'svglite', fig.ext='svg', fig.path = "SVGs/",
fig.align='center', fig.width = 7, fig.height = 3
)
```

But that looks like it's just removing svg config.

I also wasn't clear on what SETTINGS-knitr is for as it's in the
name of the file and in the call to R.

And I take it opts_chunk is from knitr and will be avaialble.

- Bob

Ben Goodrich

unread,
Jan 16, 2016, 6:24:24 PM1/16/16
to stan development mailing list
On Saturday, January 16, 2016 at 5:30:21 PM UTC-5, Bob Carpenter wrote:
I'm not sure what I'm looking at.  Which piece of it
reduces size relative to the defaults?

I think changing from .svg graphics to .png graphics is part of it, but most of it was changing the top of the .Rmd file to html_vignette instead of html_document to avoid a bunch of Bootstrapped stuff.

https://github.com/yihui/rstanarm/commit/64c71aa2c40a5ab55c8fc4c9ad27d531078e9808

Ben

Ben Goodrich

unread,
Jan 17, 2016, 8:43:30 PM1/17/16
to stan development mailing list
On Friday, January 15, 2016 at 4:26:21 PM UTC-5, Bob Carpenter wrote:
Before releasing to our users, I'd like to get some feedback
on a knitr doc I wrote as an intro to hierarchical modeling
for repeated binary trials.  It lives here:

  stan-dev/example-models/knitr/efron-morris/pool-no-pool.Rmd

Am I the only one who does not see an efron-morris directory in stan-dev/example-models/knitr?

Ben

Avraham Adler

unread,
Jan 18, 2016, 1:38:39 AM1/18/16
to stan development mailing list
Try the issue-34-case-study-egs branch.

Bob Carpenter

unread,
Jan 18, 2016, 3:30:23 PM1/18/16
to stan...@googlegroups.com
Sorry for not mentioning that. I wasn't going to
push to master until I got some review. I should've made
a pull request.

- Bob

Bob Carpenter

unread,
Jan 18, 2016, 3:40:23 PM1/18/16
to stan...@googlegroups.com
Sorry about that. I cleaned everything up and made a proper
pull request:

https://github.com/stan-dev/example-models/pull/39

- Bob

> On Jan 18, 2016, at 1:38 AM, Avraham Adler <avraha...@gmail.com> wrote:
>

Michael Betancourt

unread,
Jan 18, 2016, 8:21:54 PM1/18/16
to stan...@googlegroups.com
What needs to be installed to get rmarkdown? My copy of R doesn’t
seem to have it and nor can it find it on CRAN.

Avraham Adler

unread,
Jan 18, 2016, 8:49:23 PM1/18/16
to stan development mailing list
On Monday, January 18, 2016 at 8:21:54 PM UTC-5, Michael Betancourt wrote:
> What needs to be installed to get rmarkdown? My copy of R doesn’t
> seem to have it and nor can it find it on CRAN.
>

Does this work? <https://cran.r-project.org/web/packages/rmarkdown/index.html>

Avi

Michael Betancourt

unread,
Jan 18, 2016, 9:05:11 PM1/18/16
to stan...@googlegroups.com
This…this should not be this hard:

> rmarkdown::render("pool-no-pool.Rmd")
Error: pandoc version 1.12.3 or higher is required and was not found.
> install.packages("pandoc")
Warning message:
package ‘pandoc’ is not available (for R version 3.2.3)

Jiqiang Guo

unread,
Jan 18, 2016, 9:31:43 PM1/18/16
to stan...@googlegroups.com
Pandoc is not an r package.  See http://pandoc.org/

Michael Betancourt

unread,
Jan 19, 2016, 9:26:50 AM1/19/16
to stan...@googlegroups.com
Section: Repeated Binary Trials

I would reference n \in N in the opening sentence, otherwise the n subscripts don’t make sense until the next few sentences have been read.

“any ML conference proceedings” or “all ML conference proceedings”? :-p

Emphasize that you are using the remaining data to determine _predictive performance_.

Put “The data separates the outcome from the initial 45 at-bats from the rest of the season.” after the display of the data. Perhaps something like “The data are separated into the first 45 at bats of the season and the remaining at bats; we will use the first 45 at bats to build our model and then the remaining at bats to test the predictive performance of our model.”

Section: Pooling

There’s an awkward line break after the first two sentences.

You haven’t defined a population distribution yet so the the discussion of variances has no context here. If you wanted to include this note here, I would introduce another sentence like “Partial pooling is typically accomplished with a hierarchical or population model, where the individuals assumed to have been are drawn from a population normally distributed with mean mu and variance sigma2. When the population variance is infinite partial pooling is equivalent to no pooling, while zero population variance is equivalent to complete pooling”. Or just leave this discussion to the partial pooling section (although you don’t use a normal hierarchical model there, so the relevant of this point is questionable).

Section: Model 1: Complete Pooling

Do you want () or [] here? Typically 0 and 1 are considered valid values, in which case you want [0, 1].

“The consequences for leaving the constraint off is that the program may fail during random initialization or during an iteration because it will try values…” What is “it”? Maybe something like “…because Stan will consider mathematically invalid values outside of (0, 1)”.

The proportion sign hasn’t been used yet, why is it being described here?

“The model can be fit as follows, reading…” -> “The model can be fit as follows: <fit input>”. Note that the data are read out of the R environment by name…”

Section: Model 3: Partial Pooling (Chance of Success)

I recommend running with the default step size and adapt_delta and discuss the errors that show up. Mention how the errors indicate that Stan isn’t fitting well and we can try to tweak the step size and adapt_delta to improve the fit. It helps, but the mixing is still bad so we need a better solution (which then motivates the next section).

The Beta distribution isn’t a location-scale distribution so it is not immediately amenable to the centering/non-centering discussion. But it can be manipulated into a location-scale family by unconstraining its parameters, in which case you center/non-center the unconstrained parameters. Essentially you’ll end up with something like the model in the next section anyways.

Section: Model 4: Partial Pooling (Log Odds)

Yup, the distribution for mu and sigma will still look like a weak funnel (when the random effect variation is large, there is a weak non-identifiability between the mean and the average of the random effects). What you want to be looking at is the distribution of the random effects, alpha, verses sigma. That’s what causes all of the problems in the fit.

Section: Observed vs. Estimated Chance of Success

It would be great if more of the plot could be labeled to avoid having to describe everything in the text.

Section: Posterior Predictive Distribution

You never discuss what you can do with the predictive distribution, let alone why they provide a ranking that might approximate relative performance. This leaves a big gap when you move on to the explicit calculations.

Section: Posterior p-Values

I highly recommend a comment here about how posterior predictive p-values can’t really be calibrated (no power calculations can be done because there is no alternative hypothesis). This means that we can’t formally distinguish between a bad fit or unlikely data, and we have to be careful in how we interpret these values. Ultimately what we do in practice is assume that for small enough p-values there has to be some alternative model that fits better, and if we see such small p-values we start looking for that alternate model by expanding our null model.



On Jan 18, 2016, at 8:39 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Avraham Adler

unread,
Jan 19, 2016, 2:28:29 PM1/19/16
to stan development mailing list
IIRC, A subset of pandoc is installed automatically with rstudio, but I believe the brain trust here is primarily Linux-based, so they would need a separate install.

Avi

Bob Carpenter

unread,
Jan 19, 2016, 3:45:25 PM1/19/16
to stan...@googlegroups.com
You need to separately install pandoc on your path.

rmarkdown is not a standalone R package.

- Bob

Bob Carpenter

unread,
Jan 19, 2016, 3:47:25 PM1/19/16
to stan...@googlegroups.com
We're on a mixture of Mac and Linux. But it's
all Unix. But yes, pandoc is included as part of
Rstudio on the Mac. I had to go through what Michael
did to install it separately. Sorry --- you forget the
hassles of installs once they're in the past.

- Bob

> On Jan 19, 2016, at 2:28 PM, Avraham Adler <avraha...@gmail.com> wrote:
>
> IIRC, A subset of pandoc is installed automatically with rstudio, but I believe the brain trust here is primarily Linux-based, so they would need a separate install.
>
> Avi
>

Jiqiang Guo

unread,
Jan 19, 2016, 11:09:37 PM1/19/16
to stan...@googlegroups.com
This is long, so I have not finished.  For now, 

1) This is personal.  For the R code, I would not put semicolon at the end of statement when they are redundant. 
2) I think we cannot denote half-normal by 0.5Normal.  
3) We start with number of iteration being 10000
4) I did not see Gamma function is used, so there is no need to have 

where ...., and $\Gamma(u)$ is the [gamma


Jiqiang 

Bob Carpenter

unread,
Jan 20, 2016, 5:24:29 PM1/20/16
to stan...@googlegroups.com
1) can't help myself

2) right --- it should've been 2 * Normal

3) I think I changed the iterations to 5K from 10K,
but I'll check.

4) thanks

- Bob

Bob Carpenter

unread,
Jan 21, 2016, 5:49:31 PM1/21/16
to stan...@googlegroups.com
Thanks for catching all those inconsistencies --- some
of them came during revisions, of course. I decided to downplay
the forms of the densities, for example.

I did what you suggested everywhere; hope you don't mind that
I just included your text wholesale. So much easier than trying
to paraphrase :-)

Only two comments I couldn't follow:

> Section: Posterior Predictive Distribution
>
> You never discuss what you can do with the predictive distribution, let alone why they provide a ranking that might approximate relative performance. This leaves a big gap when you move on to the explicit calculations.

I discuss what you can do with it as soon as I define it, but
that's not near where I talk about ranking, so I don't know what
you're suggesting here.

> Section: Posterior p-Values
>
> I highly recommend a comment here about how posterior predictive p-values can’t really be calibrated (no power calculations can be done because there is no alternative hypothesis). This means that we can’t formally distinguish between a bad fit or unlikely data, and we have to be careful in how we interpret these values. Ultimately what we do in practice is assume that for small enough p-values there has to be some alternative model that fits better, and if we see such small p-values we start looking for that alternate model by expanding our null model.

Not at all sure what this means and I'm uncomfortable dropping in
a paragraph I don't understand.

Isn't it always the case for hypothesis tests that we can't distinguish
bad model vs. unlikely data? That is, we'll falsely reject null hypotheses
when we are faced with unlikely data.

- Bob

Michael Betancourt

unread,
Jan 21, 2016, 7:11:30 PM1/21/16
to stan...@googlegroups.com
>
>> Section: Posterior Predictive Distribution
>>
>> You never discuss what you can do with the predictive distribution, let alone why they provide a ranking that might approximate relative performance. This leaves a big gap when you move on to the explicit calculations.
>
> I discuss what you can do with it as soon as I define it, but
> that's not near where I talk about ranking, so I don't know what
> you're suggesting here.

So we have the sections:

Posterior Predictive Distribution
- Define the posterior predictive distribution
Evaluating Held-Out Data Predictions
- This is just a definition, there’s no motivation for what held-out data predictions can be used.
Simulating Replicated Data
- For what is this replicated data useful? Ultimately to compute posterior predictive expectations
but that’s not made clear here.
Prediction for New Trials
- Here’s where you make the jump. You talk about evaluating models, but you just assume that
held-out predictive likelihood is the right measure without any discussion as to why.

>> Section: Posterior p-Values
>>
>> I highly recommend a comment here about how posterior predictive p-values can’t really be calibrated (no power calculations can be done because there is no alternative hypothesis). This means that we can’t formally distinguish between a bad fit or unlikely data, and we have to be careful in how we interpret these values. Ultimately what we do in practice is assume that for small enough p-values there has to be some alternative model that fits better, and if we see such small p-values we start looking for that alternate model by expanding our null model.
>
> Not at all sure what this means and I'm uncomfortable dropping in
> a paragraph I don't understand.
>
> Isn't it always the case for hypothesis tests that we can't distinguish
> bad model vs. unlikely data? That is, we'll falsely reject null hypotheses
> when we are faced with unlikely data.

It is if you don’t compute a power! That’s part of the reason why significance
tests are somewhat meaningless without a power calculation. In the predictive
validation case we don’t have an alternative hypothesis so we can’t compute
anything like a power, and we’re left with the fundamental ambiguity.

Bob Carpenter

unread,
Jan 21, 2016, 9:10:31 PM1/21/16
to stan...@googlegroups.com

> On Jan 21, 2016, at 7:11 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
>>
>>> Section: Posterior Predictive Distribution
>>>
>>> You never discuss what you can do with the predictive distribution, let alone why they provide a ranking that might approximate relative performance. This leaves a big gap when you move on to the explicit calculations.
>>
>> I discuss what you can do with it as soon as I define it, but
>> that's not near where I talk about ranking, so I don't know what
>> you're suggesting here.
>
> So we have the sections:
>
> Posterior Predictive Distribution
> - Define the posterior predictive distribution
> Evaluating Held-Out Data Predictions
> - This is just a definition, there’s no motivation for what held-out data predictions can be used.

OK, I added a bit of that already in the last round
of revisions. But now I see the definition's wrong
as it was later.

It should be

log(1/M * SUM_m p(y | theta[m]))

instead of

1/M SUM_m log p(y | theta[m])

> Simulating Replicated Data
> - For what is this replicated data useful? Ultimately to compute posterior predictive expectations
> but that’s not made clear here.

I'm also plotting replications for the kind of chi-by-eye
thing Andrew likes.

> Prediction for New Trials
> - Here’s where you make the jump. You talk about evaluating models, but you just assume that
> held-out predictive likelihood is the right measure without any discussion as to why.

I wasn't trying to imply that, so I'll add some discussion
as to why it's not irrelevant (pardon me while I rev my negation
engine before the next section).


>>> Section: Posterior p-Values
>>>
>>> I highly recommend a comment here about how posterior predictive p-values can’t really be calibrated (no power calculations can be done because there is no alternative hypothesis). This means that we can’t formally distinguish between a bad fit or unlikely data, and we have to be careful in how we interpret these values. Ultimately what we do in practice is assume that for small enough p-values there has to be some alternative model that fits better, and if we see such small p-values we start looking for that alternate model by expanding our null model.
>>
>> Not at all sure what this means and I'm uncomfortable dropping in
>> a paragraph I don't understand.
>>
>> Isn't it always the case for hypothesis tests that we can't distinguish
>> bad model vs. unlikely data? That is, we'll falsely reject null hypotheses
>> when we are faced with unlikely data.
>
> It is if you don’t compute a power! That’s part of the reason why significance
> tests are somewhat meaningless without a power calculation. In the predictive
> validation case we don’t have an alternative hypothesis so we can’t compute
> anything like a power, and we’re left with the fundamental ambiguity.

I'm only vaguely familiar with power calculations from Gelman and
Hill, so I tried reading this:

https://en.wikipedia.org/wiki/Statistical_power

which says it's:

the probability that the test correctly rejects the null hypothesis (H0)
when the alternative hypothesis (H1) is true.

All the negation is like a logic puzzle.

- Bob

Michael Betancourt

unread,
Jan 21, 2016, 9:24:15 PM1/21/16
to stan...@googlegroups.com
I'm also plotting replications for the kind of chi-by-eye
thing Andrew likes.

Yeah, that’s useful (and ultimately as far as I”ll take
predictive validation) but I think it would be worth
a few sentences of explanation, like “We can simulate
data from the predictive distribution and compare it
to the measurement that we made.  If they are not
consistent then either our model is wrong or the
measurement we made is highly unlikely.  Although
we can’t formally differentiate between these two
cases, we typically assign a high prior probability
on extreme inconsistencies indicating a poor model
fit."


I'm only vaguely familiar with power calculations from Gelman and
Hill, so I tried reading this:

 https://en.wikipedia.org/wiki/Statistical_power

which says it's:

 the probability that the test correctly rejects the null hypothesis (H0)
 when the alternative hypothesis (H1) is true. 

All the negation is like a logic puzzle. 

Yeah, which is why it’s so hard to understand!
I just think of the significance as indicating
inconsistency and the power ensuring that
inconsistencies aren’t due to random variation.
That’s consistent with the more formal
definitions, but in my opinion easier to grok.

Bob Carpenter

unread,
Jan 22, 2016, 12:00:35 PM1/22/16
to stan...@googlegroups.com
Will do. And I did my homework and understand your earlier
comments now. Andrew made the point a little less
technicaly in BDA saying that it doesn't involve probabilities
of a model being true. I think I can explain that and why
we want predictive modeling and what the whole posterior
predictive thing's doing.

I made the usual mistake of understanding what was going
on myself (at least superficially in that I didn't connect
the p-values to power), but then not telling the reader.
I just went through recording what I did to reproduce all
the analyses with some definitions.

- Bob

Michael Betancourt

unread,
Jan 22, 2016, 12:03:58 PM1/22/16
to stan...@googlegroups.com
That’s why it’s so critical to start writing while you still
remember learning everything.  :-p

Jonah Gabry

unread,
Jan 22, 2016, 2:31:17 PM1/22/16
to stan development mailing list
Bob, this is awesome. I finally had a chance to read through most of it.
I think the others have made good recommendations already, but I do have one
suggestion about the R code, which is that the with() function is very useful.
There are several places where it can make your code simpler:

1) You have:

df <- data.frame(FirstName = df$FirstName,
LastName = df$LastName,
Hits = df$Hits,
At.Bats = df$At.Bats,
RemainingAt.Bats = df$RemainingAt.Bats,
RemainingHits = df$SeasonHits - df$Hits)

but you can do the same thing like this:

df <- with(df, data.frame(FirstName, LastName,
Hits, At.Bats,
RemainingAt.Bats,
RemainingHits = SeasonHits - Hits))

which avoids the df$ on every line and the repetition of variable names (which are both pet peeves of mine, but maybe that's just me).

2) The same thing goes for working with the lists that rstan::extract returns. For
example, this

phi_sim <- ss_hier$phi;
kappa_sim <- ss_hier$kappa;
df_bda3_fig_5_3 <- data.frame(x = log(phi_sim / (1 - phi_sim)),
y = log(kappa_sim));

can be done in one line without the intermediate steps of creating phi_sim and
kappa_sim:

df_bda3_fig_5_3 <- with(ss_hier, data.frame(x = log(phi / (1 - phi)),
y = log(kappa)))

3) You can do the same thing later in the file when working with the next model:

mu_sim <- ss_hier_logit$mu;
sigma_sim <- ss_hier_logit$sigma;
df_bda3_fig_5_3_logit <- data.frame(x = mu_sim,
y = log(sigma_sim));

can all be replaced with the one line

df_bda3_fig_5_3_logit <- with(ss_hier_logit, data.frame(x = mu_sim,
y = log(sigma_sim)))


Again, this really is a fantastic tutorial. And of course feel free to ignore my R code preferences. The way you have it now definitely works just fine.

Jonah Gabry

unread,
Jan 22, 2016, 2:36:53 PM1/22/16
to stan development mailing list
Oops, that last example in my previous post should use mu and sigma instead of mu_sim and sigma_sim

Bob Carpenter

unread,
Jan 24, 2016, 1:55:41 PM1/24/16
to stan...@googlegroups.com
Jonah --- how do I generalize your approach? This fails:


> df_bda3_fig_5_3_funnel <-
+ with(ss_hier_logit,
+ data.frame(x = alpha[, 1],
+ y = log(sigma)));
Error in alpha[, 1] : object of type 'closure' is not subsettableR's so confusing.


whereas this is OK:


> head(ss_hier_logit$alpha[, 1])
[1] 0.8934396 0.2556267 1.5711967 1.0863271 -0.8518395 -0.5770550


- Bob

Jonah Gabry

unread,
Jan 24, 2016, 2:52:52 PM1/24/16
to stan development mailing list
I think the approach is generalizable in the way you tried, but the error just comes from the fact that you don't actually have an "alpha" in ss_hier_logit, but rather only alpha_std (it doesn't look like alpha gets created in the Stan program). So that's why your first way doesn't work. You'd think that the second way would also fail then, but using "$" allows for partial name matching so it's accessing ss_hier_logit$alpha_std when you use ss_hier_logit$alpha. This is often convenient for interactive programming but it drives me nuts that it's the default (one of the many things about R that drives me nuts).

Jonah

Bob Carpenter

unread,
Jan 24, 2016, 3:11:38 PM1/24/16
to stan...@googlegroups.com
Thanks. It worked when I named the variable properly.

Personally, I'd prefer completion to auto-extension in an
interactive setting.

- Bob

> On Jan 24, 2016, at 2:52 PM, Jonah Gabry <jga...@gmail.com> wrote:
>
> I think the approach is generalizable in the way you tried, but the error just comes from the fact that you don't actually have an "alpha" in ss_hier_logit, but rather only alpha_std (it doesn't look like alpha gets created in the Stan program). So that's why your first way doesn't work. You'd think that the second way would also fail then, but using "$" allows for partial name matching so it's accessing ss_hier_logit$alpha_std when you use ss_hier_logit$alpha. This is often convenient for interactive programming but it drives me nuts that it's the default (one of the many things about R that drives me nuts).
>
> Jonah
>

Jonah Sol Gabry

unread,
Jan 24, 2016, 3:33:41 PM1/24/16
to stan...@googlegroups.com
Yeah, and RStudio (maybe R in general, I forget) will bring up suggestions if you type x$partial and there's an element in x whose name starts with "partial". But instead of making you use one of the suggestions (or letting it autocomplete) you are allowed to just use x$partial. This makes it quite dangerous to use "$" when writing code for an R package (which I've learned the hard way a few times). 

There are also some weird rules about when the names are too ambiguous for the partial matching, but it's hard to remember those and the documentation is pretty abysmal.

You received this message because you are subscribed to a topic in the Google Groups "stan development mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-dev/Cea1RdQZMSA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-dev+u...@googlegroups.com.

Krzysztof Sakrejda

unread,
Jan 24, 2016, 3:45:06 PM1/24/16
to stan development mailing list
On Sunday, January 24, 2016 at 3:33:41 PM UTC-5, Jonah Gabry wrote:
> Yeah, and RStudio (maybe R in general, I forget) will bring up suggestions if you type x$partial and there's an element in x whose name starts with "partial". But instead of making you use one of the suggestions (or letting it autocomplete) you are allowed to just use x$partial. This makes it quite dangerous to use "$" when writing code for an R package (which I've learned the hard way a few times). 

It's nothing really with RStudio, just R in general! It's another one of those charming quirks... except R will also do partial matching on argument names in functions (!) so there's no way to get away from it.

Krzysztof

Krzysztof Sakrejda

unread,
Jan 26, 2016, 9:03:41 AM1/26/16
to stan development mailing list
On Friday, January 15, 2016 at 4:26:21 PM UTC-5, Bob Carpenter wrote:
> Before releasing to our users, I'd like to get some feedback
> on a knitr doc I wrote as an intro to hierarchical modeling
> for repeated binary trials. It lives here:
>
> stan-dev/example-models/knitr/efron-morris/pool-no-pool.Rmd
>
> As usual, I had a hard time judging what to explain and what
> to take for granted. I tried to focus on Stan. Any tips on
> making the RStan better would be appreciated. Andrew, now's the
> time to jump in on ggplot2 defaults! Or I can ask you in person
> and we can create a theme you can live with.
>
> I'm happy to send the 2.8MB HTML file it produced (encapsulated)
> if you'd rather not run it yourself. (If you do want to run it,
> you only need:
>
> > rmarkdown::render("pool-no-pool.Rmd")

I got as far as this step (just did it yesterday) and go this:

label: unnamed-chunk-14 (with options)
List of 1
$ comment: logi NA

|........................ | 37%
ordinary text without R code

|......................... | 38%
label: unnamed-chunk-15
Quitting from lines 801-813 (pool-no-pool.Rmd)
Error in with(ss_hier_logit, data.frame(x = alpha_1_centered, y = log(sigma))) :
object 'ss_hier_logit' not found


>
> to get the basic HTML. I couldn't get it to keep the figures
> separate, even trying what Ben suggested (I'm probably doing it
> wrong).
>
> Here's the abstract:
>
> This note illustrates the effects on posterior inference of pooling
> data (aka sharing strength) across items for repeated binary trial
> data. It provides Stan models and R code to fit and check predictive
> models for three situations: (a) complete pooling, which assumes each
> item is the same, (b) no pooling, which assumes the items are
> unrelated, and (c) partial pooling, where the similarity among the
> items is estimated. We consider two hierarchical models to estimate
> the partial pooling, one with a beta prior on chance of success and
> another with a normal prior on the log odds of success. The note
> explains with working examples how to (i) fit models in RStan and plot
> the results in R ussing ggplot2, (ii) estimate event probabilities,
> (iii) evaluate posterior predictive densities to compare model
> predictions, (iv) rank items by chance of success, (v) perform
> multiple comparisons in several settings, (vi) replicate new data for
> posterior p-values, and (vii) perform graphical posterior predictive
> checks.
>
> - Bob

Bob Carpenter

unread,
Jan 26, 2016, 5:06:42 PM1/26/16
to stan...@googlegroups.com
OK, it's all revised and checked in if anyone wants
a final gander. Thanks for all the help --- it's very
nice to get world-class feedback from multiple people on my
homework!

I didn't go to the next level of trying to write a dissertation
on PPCs as part of an intro tutorial.

- Bob
Reply all
Reply to author
Forward
0 new messages