improving doc for posterior_linpred with new levels

[wds15]

It seems that posterior_linpred samples new levels of random effects only once per draw even if the new level varies and as such should be re-drawn multiple times.

Here is a minimal example:

library(brms)
library(mvtnorm)
library(dplyr)

set.seed(425435)

R <- matrix(c(1, 0.1, 0.1, 1), 2, 2, byrow=TRUE)
R

J <- 20
group_draw <- matrix(rmvnorm(J, sigma=R), ncol=2, dimnames=list(NULL, list("a", "b"))) %>%
    as.data.frame() %>%
    mutate(id=1:J)

group_draw

fake  <- expand.grid(id=1:J, x=seq(0,1,length=11)) %>%
    left_join(group_draw) %>%
    arrange(id, x) %>%
    mutate(y_hat= a + b * x, y = rnorm(n(), y_hat), id=factor(id))

head(fake)

fake

form <- bf(y ~ 1 + x + (1+x|id), center=FALSE)

get_prior(form, data=fake, family=gaussian)

model_prior <- prior(normal(0, 2), class=b, coef=Intercept) +
    prior(normal(0, 1), class=b) +
    prior(normal(0, 1), class=sd, coef=Intercept, group=id) +
    prior(normal(0, 0.5), class=sd, coef=x, group=id)

fit  <- brm(form, data=fake, prior=model_prior, family=gaussian, cores=4)


fit

df  <- data.frame(id=factor((J+1) : (2*J)), x=c(0))

df

plin1  <- posterior_linpred(fit, newdata=df, allow_new_levels=TRUE)
plin2  <- posterior_linpred(fit, newdata=df, allow_new_levels=TRUE)


plin1[1,]
plin2[1,]

The output is:

> plin1[1,]
 [1] 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684
 [8] 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684
[15] 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684 0.1073684
> plin2[2,]
 [1] -0.1205381 -0.1205381 -0.1205381 -0.1205381 -0.1205381 -0.1205381
 [7] -0.1205381 -0.1205381 -0.1205381 -0.1205381 -0.1205381 -0.1205381
[13] -0.1205381 -0.1205381 -0.1205381 -0.1205381 -0.1205381 -0.1205381
[19] -0.1205381 -0.1205381

So the random effect which varies by id is kept constant for an entire draw, but that's just not right as the id varies across the data set which I wish to simulate.

Is this a user error or a bug in brms?

The above code ran with brms 2.13.0, but it is the same misery in 2.10.0

As the new data set as many new id's I would expect that I get a newly drawn random effect for each row of my data set.

[wds15]

rstanarm seems to behave the same way here.

[paul-buerkner]

This has to do with the sample_new_levels argument. Does it work as you expect, when you set sample_new_levels = "gaussian"?

[wds15]

Using sample_new_levels = "gaussian" makes it look a lot better. So it looks like what I want.

Would it be possible to add a hint in the documentation of posterior_predict to the allow_new_levels option that behind the "..." given to extract_draws there is this important option?

Thanks for the quick response and sorry for the noise... but I was really confused from this; glad it's just another switch...

[wds15]

Let's change the issue to what this should be. Improving the doc would be great. So wherever allow_new_levels pops up, a hint to the sample_new_levels should be made, I think. That would maybe have saved my afternoon yesterday and hence possibly others as well.

[paul-buerkner]

But isn't sample_new_levels documented in the same place as allow_new_levels? Can you clarify where you would like to see the doc updated?

[wds15]

I think all the predict functions should probably get a new section about "Predicting new levels" which you can include with roxygen template stuff.

The posterior_linpred function does not even mention allow_new_levels nor does it mention sample_new_levels. The posterior_predict also does not mention these super important parameters at all. They are hidden behind the ....

EDIT: The docs do say something about prepare prediction functions, but that is not clear enough what that all entails - at least that was not clear for me.

[paul-buerkner]

Yeah, you are right those arguments are hard to find. I added a paragraph to the details section pointing to these arguments.

[wds15]

Thanks for fixing the doc.

Just wanted to note that I still find the current default way super-dangerous. I mean, what I did is to simulate new random effect levels which do vary from row to row in my data-set... but what I get is just one realisation of the random effect per draw (thus being constant across the entire data-set). I don't understand the default right now, but I am really not sure if users would ever want that. There should be at least a warning being emitted as the model is structurally being wrongly simulated to my understanding.

[paul-buerkner]

I understand your concern. Since rstanarm does mostly the same thing here, I would like to hear the opinion of @jgabry and @bgoodri on this issue as well.

[bgoodri]

With rstanarm, we have to put an extra level into the MCMC part (it gets multiplied by zero) in order to get draws from its posterior distribution. We might be able to make that more sane with a standalone generated quantities block. In any event, brms shouldn't feel constrained by that.