Why am I getting noisy predictions from brms::predict?


First, thank you for BRMS and all the underlying work. I’m hoping to get some clarification of why my predictions are noisy when I would expect them to be linear? (sorry if I’m missing something obvious) I am working on a multi-output model but my questions is focused on one of the levels which would apply to all ultimately. First, my data is a time series data set of health spending from many countries with many years of missing data in some countries and nearly complete data in others. My base model simply fits the log of spending with year.

`ln(health_exp) ~ year

fit <- brm(ln_value_pc ~ year, data = data_small_total)

Family: gaussian
Links: mu = identity; sigma = identity
Formula: ln_che_pc ~ year
Data: data_small_total (Number of observations: 1340)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
ICs: LOO = NA; WAIC = NA; R2 = NA

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 29.90 18.16 -5.75 64.87 4000 1.00
year -0.01 0.01 -0.03 0.01 4000 1.00

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 1.88 0.04 1.82 1.96 3747 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

I would expect to simply predict a linear trend when I use predict to fill in all the years of data across my years. However, when using predict or posterior.predict the within country annual variation does not follow a simply linear trend. Below is the head of these predictions:

pred_che_fit <- predict(fit_che, newdata = data_small_total, allow_new_levels = TRUE, sort = TRUE)

Estimate Est.Error 2.5%ile 97.5%ile
1: 7.306073 1.411603 4.515836 9.929672
2: 7.343608 1.411364 4.573684 10.066197
3: 7.345000 1.411676 4.510049 10.068934
4: 7.359616 1.415086 4.641277 10.099347
5: 7.323695 1.410809 4.547681 10.137092
6: 7.349641 1.409937 4.639661 10.217046
Again, my interest is why are the above estimates and percentiles going up and down over year when it should be linear?

Thanks in advance! Matt

  • Operating System: Using Rstudio on a multi core cluster computer
  • brms Version: 2.2.0

A posterior prediction, by definition, adds noise to the conditional mean function. If your interest lies in how a predictor influences the conditional mean function, then you could investigate the output of posterior_linpred (with transform = TRUE in GLMs). But most people that come from a frequentist or supervised learning background confuse the latter with the former.

1 Like

Thank you Ben. This does make the predictions linear. As a Bayesian novice, are there any good background documents that would explain the underlying differences between posterior_linpred and posterior_predict? Thanks again, Matt

In the case of Gaussian outcomes, posterior_predict is equal to posterior_linpred plus noise with standard deviation given by the posterior distribution of \sigma. So, posterior_predict is what your model thinks the distribution of the outcome should look like for given values of the predictors. Or if you are just evaluating it at the observed values of the predictors, then it is supposed to look like the marginal distribution of the outcome.

hi !

posterior_linpred(fit, transform = TRUE)
seem to both give
= \{E[y_i | x_i, \theta^{(s)}]\}_{s,i}
An S x N matrix of draws from the expected outcome, where S is the number of posterior draws, N is the number of observations.

Taking the mean across posterior draws for each column (observation) gives:
= \Big \{E \Big [E[y_i | x_i, \theta^{(s)}]\ \Big |\ x_i, \text{data} \Big ] \Big \}_{i} = \{E [y_i | x_i, \text{data} ] \}_{i}
the means of the posterior predictive distribution.

The documentation for posterior_epred(fit) is a bit confusing because it says the output are these means themselves, as if the averaging over posterior draws were already done.

Does the above seem right ?

Thanks so much !