rstanarm: how to specify a varying slope and intercept hierarchical model

[Luke Wiklendt]

I've created a simple synthetic hierarchical dataset, with 5 groups of 2D
points, and a linear relationship between variables x and y.
The true intercepts for each group are [6, 8, 10, 12, 14] and all slopes
are -1, and some noise added (see circles in plots).
The titles of the figures represent the model passed to stan_glmer. The
first "y ~ x + (1|group)" seems to work very well, however, I can't make
sense of what is happening to the "y ~ x + (x|group)" model.
If I plot the overall samples of intercept and slope (bottom two plots in
grey), they seem to work fine for both, it's just the group intercepts and
slopes in the "y ~ x + (x|group)" that I don't understand.

[Andrew Gelman]

1. Those are pretty graphs. You might also try an array of 5 little plots next to each other, in case all the color overlap gets confusing.
2. Results of varying-intercept, varying-slope models will depend on how x is centered. Is a 0 value of x reasaonable in your example?
A

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

[Ben Goodrich]

On Wednesday, February 3, 2016 at 7:16:47 PM UTC-5, Luke wrote:

If I plot the overall samples of intercept and slope (bottom two plots in
grey), they seem to work fine for both, it's just the group intercepts and
slopes in the "y ~ x + (x|group)" that I don't understand.

Just to make sure we are on the same page here, I believe it is coef() function that returns what is of interest to you, which is essentially the sum of the output of fixef() and ranef() by group level. Is that what is being plotted?

Ben

[Lukasz Wiklendt]

That might be the issue. I'm not using the coef() function, I'm extracting the samples of the variables and combining the values myself.

I'm not very familiar with r, I'm running rstanarm from python via rpy2. I obtain the samples by running

as.matrix(fit$stanfit)

and the variable names with

names(fit$stanfit)

where the fit object was obtained by the call to sample

fit = stan_glmer(...)

I zip up the names and matrix and pass it to a pandas dataframe, lets call it "df".

What I've plotted given a sample index i, an overall (grey) line specified by something to the effect of

df['(Intercept)'][i] + df['x'][i] * np.array([-1, 7])


and a per-group (coloured), say for group “a", a line specified by

df['(Intercept)'][i] + df['b[(Intercept) group:a]'][i] + (df['x'][i] + df['b[x group:a]'][i]) * np.array([-1, 7])


whereas for the varying-intercept-only model each per-group line is

df['(Intercept)'][i] + df['b[(Intercept) group:a]'][i] + df['x'][i] * np.array([-1, 7])


given the following summary output (showing the names of the sampled variables)

                                5%        95%
(Intercept)             8.52958750 11.2192503
x                      -1.07276581 -0.8569989
b[(Intercept) group:a] -5.28344691 -2.5955489
b[(Intercept) group:b] -0.14002830  0.1422533
b[(Intercept) group:c] -3.39296111 -0.7066493
b[(Intercept) group:d] -0.07888252  0.1730402
b[(Intercept) group:e] -1.31888695  1.3628254
b[x group:a]            0.70504964  3.4770025
b[x group:b]           -0.14329447  0.1276788
b[x group:c]            2.49437147  5.2757959
b[x group:d]           -0.11953980  0.1450246
b[x group:e]           -2.79710230  2.8382719
sigma                   0.10347701  0.1544424


Luke

--
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/E-7Zau2-Jno/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Thursday, February 4, 2016 at 12:24:36 AM UTC-5, Luke wrote:

That might be the issue. I'm not using the coef() function, I'm extracting the samples of the variables and combining the values myself.

I would use coef(). If you are looking at b[(Intercept) group:?] and b[x group:?] those represent the deviation of group ? from the common intercept or slope respectively. The stan_[g]lmer() functions in rstanarm use the same parameterization that is used in the [g]lmer() functions in lme4, which is described in section 1.1 of a recent article by the lme4 authors

http://www.jstatsoft.org/article/view/v067i01

Ben

[Luke]

There still seems to be a problem. I've created a cut-and-paste script in R to reproduce the problem, excuse it's crudity, it's the first script in R I've written.
The model from y ~ x + (1|grp) works fine, but the model y ~ x + (x|grp) shows a problem visible by the histogram of posterior marginals and from the result of the coef() function.
Running the model y ~ x + (x|grp) with lme4's lmer function and sending the result to coef() shows a correct fit.

library(rstanarm)
library(ggplot2)
library(lme4)
set.seed(0)
num_grps <- 5
num_pnts_per_grp <- 10
x <- matrix(rnorm(num_pnts_per_grp * num_grps, 1:num_grps, 0.5), nrow=num_grps)
y <- 1:num_grps * 2 + 4 + rep(-1, num_grps) * x + rnorm(x, sd=0.1)
grp <- rep(letters[1:num_grps], num_pnts_per_grp)
df <- data.frame(grp, x=c(x), y=c(y))
p <- ggplot(df, aes(x, y))
p + geom_point(aes(color = grp))

fit <- stan_lmer(y ~ x + (x|grp), data=df, seed=0, chains=8, cores=8, iter=1000)
lmfit <- lmer(y ~ x + (x|grp), data=df)

plot(fit, plotfun="hist", pars=names(fit$stanfit))

coef(fit)
coef(lmfit)

Luke

[Jonah Gabry]

Thanks Luke, that script is helpful. I'll hopefully have time to dig into this today and see what's going on.

Jonah

[Jonah Gabry]

Hi Luke, 

If I run the models in your script and then call ranef() I get 

ranef(fit)

$grp
  (Intercept)     x
a       -3.88  0.00
b       -1.88 -0.03
c       -0.24  1.76
d        0.03  4.31
e       -0.08  0.04

and 

ranef(lmfit)

$grp
  (Intercept)     x
a       -3.94  0.01
b       -1.93 -0.03
c       -0.30  0.10
d        1.76  0.03
e        4.41 -0.11

I think what's happening is that you're getting the right estimates but some indices are getting screwed up when reporting them. For example, if you look at the x estimates for groups c and d in the first set of results and the (Intercept) estimates for groups d and e in the second set of results they match. Looking directly at the posterior draws also leads me to believe that the correct model is being estimated and the issue is just with labeling. We'll get this sorted out as soon as possible. Thanks for bringing this to our attention!

[Paul Buerkner]

I stumbled upon the same (or a related) problem a few days ago. In my example, the first half of the random effects estimates (intercept and slope) were correct and the second half was different than the results of lme4 and brms. Also WAIC (and LOO) estimates were incorrect I believe. If you need another reproducable example for this problem, I can post this here or on github.

[Jonah Sol Gabry]

Another reproducible example would be helpful. Thanks Paul.

--

[Bob Carpenter]

If you do full Bayes on these models and use posterior means or
medians to estimate, you wouldn't expect the same result as lme4
would you? Or did you guys implement MML in RStanARM?

Is there a writeup somewhere of the expression language for models
and the algorithms you're using?

- Bob

> --

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

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

[Jonah Sol Gabry]

On Thu, Feb 4, 2016 at 2:36 PM, Bob Carpenter <ca...@alias-i.com> wrote:

If you do full Bayes on these models and use posterior means or
medians to estimate, you wouldn't expect the same result as lme4
would you?  Or did you guys implement MML in RStanARM?

For lme4-style models we do HMC or either of the VB algorithms, not MML. It's true that we don't expect the same results, but in this case there's definitely something getting mixed up when we're renaming parameters (to give them more informative names than they have in the Stan code) before returning the results to the user. I've verified that it's only happening with a subset of lme4-style models. And for those models the estimates from Stan look good, there are a just a few misnamed parameters from the user's perspective.

 

Is there a writeup somewhere of the expression language for models
and the algorithms you're using?

Do you mean for the lme4-style language?

[Paul Buerkner]

Here we go :-)

[Ben Goodrich]

On Thursday, February 4, 2016 at 2:37:53 PM UTC-5, Bob Carpenter wrote:

If you do full Bayes on these models and use posterior means or
medians to estimate, you wouldn't expect the same result as lme4
would you?  

Strictly speaking, no, but the default prior on the group-specific parameters is pretty weak. In this case, we did get something similar but messed up the indices in the output.
 

Or did you guys implement MML in RStanARM?

No
 

Is there a writeup somewhere of the expression language for models
and the algorithms you're using?

Not really, we're mostly just relying on foo's syntax in stan_foo. The lme4::glmer and gamm4::gamm4 syntaxes are the only ones that are tricky

http://www.jstatsoft.org/article/view/v067i01

has the syntax for glmer and I myself don't really understand how gamm4 extends the spline stuff onto the lme4 syntax.

The possible algorithms are NUTS (always), meanfield / fullrank ADVI (except for stan_lm where it didn't work well in tests), and LBFGS (except for stan_lm, stan_glmer --- because plain optimization doesn't make sense --- and things that call stan_glmer like stan_lmer, stan_gamm4, stan_glmer.nb, etc.). Those are briefly explained in help("rstanarm-package").

Ben

[Bob Carpenter]

Thanks --- that answers both of my questions.

- Bob

[Jonah Gabry]

Hi Luke, 

This should be fixed now in the development version on GitHub (the readme has installation instructions). The fix will also be in the next release that we're planning to get out this week. Thanks again for reporting the problem. 

Jonah

[John Kastellec]

Hi, 

I seem to be  having a similar problem as Luke. I tried re-installing rstanarm using the development version, but got the following error.

Thanks,

John

> install_github("stan-dev/rstanarm", args = "--preclean", build_vignettes = TRUE)

Downloading GitHub repo stan-dev/rstanarm@master

Installing rstanarm

'/Library/Frameworks/R.framework/Resources/bin/R' --no-site-file --no-environ --no-save --no-restore CMD build  \

  '/private/var/folders/0z/jxhfxg4j1dnf1tzd8y7d3pk83y7p1v/T/RtmpbNAg7L/devtools97591311b4/stan-dev-rstanarm-d2058cf'  \

  --no-resave-data --no-manual 

* checking for file ‘/private/var/folders/0z/jxhfxg4j1dnf1tzd8y7d3pk83y7p1v/T/RtmpbNAg7L/devtools97591311b4/stan-dev-rstanarm-d2058cf/DESCRIPTION’ ... OK

* preparing ‘rstanarm’:

* checking DESCRIPTION meta-information ... OK

* cleaning src

* running ‘cleanup’

* installing the package to build vignettes

* creating vignettes ... ERROR

Warning in engine$weave(file, quiet = quiet, encoding = enc) :

  Pandoc (>= 1.12.3) and/or pandoc-citeproc is not available. Please install both.

Error: processing vignette 'aov.Rmd' failed with diagnostics:

It seems you should call rmarkdown::render() instead of knitr::knit2html() because aov.Rmd appears to be an R Markdown v2 document.

Execution halted

Error: Command failed (1)

On Wednesday, February 3, 2016 at 7:16:47 PM UTC-5, Luke wrote:

[Jonah Gabry]

Hi John, 

We're working to get an update on cran ASAP, but in the meantime if you want to install from GitHub I think you can get around the issue you're having in one of three ways. It looks like you don't have pandoc installed, which is needed to build the vignettes, so you can either (1) install pandoc yourself, (2) use RStudio (which I believe comes with pandoc), or (3) install rstanarm from GitHub without the build_vignettes argument (or setting it to FALSE).

Hope that helps, 

Jonah 

[John Kastellec]

Hi,

Option (3) worked--thanks.

But I still think the labels on the random effects might be off. Specifically, using the a stan_glmer model, I first obtained the predicted values from fitted(model). I then calculated the predicted values by hand.  The two differ substantially.  

I ran the same model using plain glmer and a) the predicted values from using predict() correlate at .99 across the glmer and rstanarm stan model, so I think the predict() function from the latter is working properly. I also was able to recover the predicted values manually using the results of the glmer model.

I could post a replicable example, if it helps.

Thanks,

John

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

[Jonah Sol Gabry]

Hi John,

Yeah if you could share a reproducible example that would be very helpful. Thank you!

Jonah

[John Kastellec]

Hi,

Please see the zip file here: https://dl.dropboxusercontent.com/u/8802928/kastellec_replicate.zip

Thanks very much.

Best,

John

[Ben Goodrich]

On Saturday, February 13, 2016 at 4:30:50 PM UTC-5, John Kastellec wrote:

Please see the zip file here: https://dl.dropboxusercontent.com/u/8802928/kastellec_replicate.zip

I think you may have installed a version where this was only half-fixed. When I do it with the current version on GitHub, I get a correlation of 1 between fitted() and invlogit(...) for both lme4 and rstanarm. Can you try it again?

Thanks,
Ben

[John Kastellec]

Hi Ben,

I tried to update, but am still having the same problem.  I have version 2.9.0-3 installed -- is that the latest?

Thanks,

John

[Ben Goodrich]

On Monday, February 15, 2016 at 10:49:09 AM UTC-5, John Kastellec wrote:

I tried to update, but am still having the same problem.  I have version 2.9.0-3 installed -- is that the latest?

There are many commits on GitHub where the version number is 2.9.0-3 (including the one now on most CRAN mirrors). Another issue is that you have to make sure that your workspace does not have any objects created by an older version of rstan[arm] when you update and / or when you load the rstanarm package.

Ben

[Jonah Gabry]

Hi John, thanks a lot for the example. Sorry for the delayed response, I was quite ill over the weekend. Regarding your results, I'm seeing the same thing as Ben (correlation of 1), so I think Ben's right that it must be an issue with the installation or the wrong version being loaded, etc. It looks like the binaries have built for the new CRAN release so you could try installing from CRAN (maybe set repos="https://cran.rstudio.com/" in your call to install.packages in case you have a different default mirror where the binaries aren't built yet). But, like Ben suggested, definitely also double check if your workspace has anything from a previous version.

Jonah

[John Kastellec]

Hi,

Thanks for the quick responses. I tried your suggestions and am still getting the error. The STAN model object itself is based on an old versions of rstanarm. But I assume that can't be the issue if you are getting the correct results?

Thanks,

John

On Mon, Feb 15, 2016 at 12:12 PM, Jonah Gabry <jga...@gmail.com> wrote:

Hi John, thanks a lot for the example. Sorry for the delayed response, I was quite ill over the weekend. Regarding your results, I'm seeing the same thing as Ben (correlation of 1), so I think Ben's right that it must be an issue with the installation or the wrong version being loaded, etc. It looks like the binaries have built for the new CRAN release so you could try installing from CRAN (maybe set repos="https://cran.rstudio.com/" in your call to install.packages in case you have a different default mirror where the binaries aren't built yet). But, like Ben suggested, definitely also double check if your workspace has anything from a previous version.

Jonah

--

[Jonah Sol Gabry]

Hi John, when you say the Stan model object is based on old versions what do you mean? Do you the fitted model object or the underlying compiled Stan model that's hidden from the user? When I checked the results I refit the models using the new rstanarm so there wasn't anything from a previous version. But I might be misunderstanding you.

Jonah

[John Kastellec]

Ah, I misunderstood. I thought the issue was in the post-processing, so I had not run the model again. I just did and everything works.  Thanks again.

Best,

John

[Jonah Sol Gabry]

Ok great, sorry for the confusion. Glad it's working finally!

[Jonah Sol Gabry]

And just to clarify, because I think I see the confusion now, it's still true that the model was being estimated correctly in the previous version. But the labelling mistake was occurring immediately after estimation when renaming the parameters (changing their names from those used in the underlying Stan code to the more substantive names returned to the user). So it was happening after estimation but before returning the fitted model object. That is, the fitted model object had the right estimates in it, they were just labeled incorrectly. The post-estimation functions were working properly, but they assumed the labels were correct. 

Anyway, I'm glad we got that sorted out. Thanks for your patience!

Jonah

[Ben Goodrich]

On Tuesday, February 16, 2016 at 3:44:36 PM UTC-5, Jonah Gabry wrote:

And just to clarify, because I think I see the confusion now, it's still true that the model was being estimated correctly in the previous version. But the labelling mistake was occurring immediately after estimation when renaming the parameters (changing their names from those used in the underlying Stan code to the more substantive names returned to the user). So it was happening after estimation but before returning the fitted model object. That is, the fitted model object had the right estimates in it, they were just labeled incorrectly. The post-estimation functions were working properly, but they assumed the labels were correct. 

Anyway, I'm glad we got that sorted out. Thanks for your patience!

And thus everyone needs to install the new version of rstanarm from CRAN and redo any lme4-style model from a previous version of rstanarm (stan_lmer, stan_glmer, stan_glmer.nb, stan_gamm4) that has group-specific coefficients.

Ben