R linear model interface
Richard McElreath
On Wed, Nov 12, 2014 at 10:06 AM, Bob Carpenter wrote:
[added Ben and Jiqiang and Andrew to cc list --- they're
our R experts and lmer experts]
Richard --- are you on the dev list? I'd rather have the
discussion there if you don't mind.
This is interesting. I had guessed that lm/lmer, glm/glmer, would just
use R's built-in parser, since the description on the R builtins seem
to have all the syntax.
There are two paths I could see us taking, because of license
issues:
1. Use pre-existing GPL-licensed components (from R, lme, etc.)
and then drop those into RStan (but not the other interfaces).
2. Write our own parser, BSD license it, and drop it into Stan so
it'll be available from all of our interfaces. It would probably take
me two weeks or so to write the parser to take the language to an
abstract syntax tree.
One benefit of approach 2 is that Stan's variables are a proper
subset of R's (we don't allow "." in identifiers). On the other
hand, we could probably take approach 2 and just replace "." with
"_" and hope there are no conflicts (seems unlikely, but it'd be
nice to flag).
I'm not sure what the tradeoffs are of creating a large data matrix.
I think it's going to be wasteful computationally in both speed
and memory. But maybe it would be worth it for the clarity?
- Bob
If we want to make this RSt
> On Nov 12, 2014, at 9:31 AM, Richard McElreath wrote:
>
> Hi, guys.
>
> I'd be thrilled to help out. For your question about formal specifications, I haven't seen any. Some of the conventions are...weird. For example, subtracting the intercept with -1 or +0. My parser handles most things correctly, provided they are in conventional order, but I do have some problem cases flagged (one or two on github). I don't think supporting everything is really necessary, though.
>
> I see a few routes to making your path easier.
>
> (1) glmer2stan has a linear formula parser. It needs to be updated to handle more "modern" Stan code, but it is functional. That said, it isn't a great compilation algorithm, so would be awful to maintain going forward. I hacked it together as an experiment, really. As usual with such coding experiments, by the end I knew a better approach. But the parser is definitely a good start.
>
> (2) My 'rethinking' package has another formula parser (the same one, basically) and another model compilation algorithm. It uses an intermediate model specification layer, the 'map' formula lists that are in my in-progress book. So the function glimmer() takes lmer-style formulas and makes more math-stats style lists of distributional assumptions. e.g.:
>
> > library(rethinking)
> > library(lme4)
> > data(sleepstudy)
> > glimmer( Reaction ~ Days + (Days | Subject), sleepstudy )
> alist(
> Reaction ~ dnorm( mu , sigma ),
> mu <- Intercept +
> b_Days*Days +
> v_Intercept[Subject] +
> v_Days[Subject]*Days,
> Intercept ~ dnorm(0,10),
> b_Days ~ dnorm(0,10),
> c(v_Intercept,v_Days)[Subject] ~ dmvnorm2(0,sigma_Subject,Rho_Subject),
> sigma_Subject ~ dcauchy(0,2),
> Rho_Subject ~ dlkjcorr(2),
> sigma ~ dcauchy(0,2)
> )
>
> Then those formula lists can be compiled to Stan code using map2stan(). map2stan is templated, and so very expandable. I plan to maintain rethinking going forward, so may be easiest to just pull out the glimmer/map2stan functions. map2stan needs a "Cholesky turbo" option, and it needs to detect your embedded link function distributions (poisson_log_log e.g.), but otherwise I believe it's up-to-date with RStan 2.5.
>
> (3) It might make more sense for your purposes to just hijack lme4's parser and then use it's internal matrix representation of the models to build the Stan code. That way, you get better insurance of compatibility with what users expect from lme4, and you get all the predictors into matrix form right away. lme4 is GLPv3 (like my code), so no obstacles there, I don't believe. I actually used bits of lme4's hidden methods to get my parser working, but I parsed into a different representation, so results are not always identical.
>
> So in sum, thrilled to help out as best I can. Let me find time soon (probably Thursday at earliest) to send you a short script with examples of what rethinking's glimmer|map2stan can do right now.
>
> You can start poking around my parser, from the github code. Check out:
>
> https://github.com/rmcelreath/glmer2stan/blob/master/R/parse_formula.R
>
> or
>
> https://github.com/rmcelreath/rethinking/blob/master/R/glimmer.R
>
> And yes, all my stuff should be GPLv3.
>
> Best,
> Richard
>
> On Tue, Nov 11, 2014 at 10:48 AM, Rob Trangucci wrote:
> Hi Richard,
>
> I recently joined the Stan team, and I'm working on multilevel models for survey data. Many potential Stan users are more familiar with lmer modeling syntax rather than Stan, so implementing an "lmer2stan" function in Rstan is one of my priorities.
>
> After some digging, I see that you've already done this (and more) with your glmer2stan package! If possible, I'd like to extend what you've already done by optimizing the Stan code glmer2stan generates to conform to our current understanding of best practices for multilevel models, and perhaps bring it under the RStan umbrella.
>
> Let me know if this is OK with you, and if so we can talk more. One of our central questions is whether there are any formal language specifications for the modeling language in lmer (the closest I could find for R's modeling language is in Wilkinson and Rogers' "Symbolic Description of Factorial Models for Analysis of Variance").
>
> Best,
> Rob
>
Richard McElreath
My parser does the same thing, but uses it's own parsed representation. You'll see in my github code that there are internal function nobars() and findbars() that recursively crawl through the tree structure of the formulas (R formulas are really function call trees) to do the splitting. Those are lme4 functions, for the most part. As much as I enjoy recursion, in hindsight it would probably be easier to convert the whole formula to text and use grep to do the splitting.
Bob's suggestion #2 sounds ambitious and very useful. Would be fantastic for CmdStan to have a linear formula parser.
Bob Carpenter
> On Nov 13, 2014, at 9:20 AM, Richard McElreath <rmcel...@gmail.com> wrote:
>
> What I discovered when I started glmer2stan was that lme4's use of the parentheses for varying effects is idiosyncratic. Other and older packages, like lme and MCMCglmm, use separate formulas for fixed and varying effects. So lme4 does its own parsing into separate formulas, and then it passes formulas to the internal methods model.matrix (and model.frame I think) for processing. But it has to split apart the formulas first.
>
> My parser does the same thing, but uses it's own parsed representation. You'll see in my github code that there are internal function nobars() and findbars() that recursively crawl through the tree structure of the formulas (R formulas are really function call trees) to do the splitting. Those are lme4 functions, for the most part. As much as I enjoy recursion, in hindsight it would probably be easier to convert the whole formula to text and use grep to do the splitting.
>
> Bob's suggestion #2 sounds ambitious and very useful. Would be fantastic for CmdStan to have a linear formula parser.
of our doing something reasonable for RStan in the short term.
Thanks much for the rundown on lme4 parsing. I wonder
if JAGS does the same thing?
One reason we were considering writing our own parser is that
might want to extend the language somehow, particularly to allow
the representation of the models in Gelman and Hill's book which
can't be expressed in the lme4 language.
As to matrices, consider something like the simplest model. It's
basically an (M + N) dimensional regression where there are M
students and N questions. But only two of the dimensions are
ever on for any example, so coding this using matrix multiplication
would be really costly in Stan (we don't have a good way to optimize
out all the zeros, so they all get autodiffed).
As for extensions, if you consider the standard 2PL version of IRT,
you get another N parameters, but it's no longer a simple linear
model, so I think you're beyond what lme4 can do. There's also
the identifiability issue.
- Bob