# Predicting the variance, taking the mean into account

**URL:** <https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955>\
**Category:** brms\
**Tags:** posterior-predictive, brms\
**Created:** [October 17, 2024, 5:29pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955 "2024-10-17T17:29:30Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ax3man](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ax3man/32/11894_2.png) [@Ax3man](https://discourse.mc-stan.org/u/Ax3man)\
**Post date:** [October 17, 2024, 5:29pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/1 "2024-10-17T17:29:30Z")

</div>

I am modelling the variance among individuals for a trait by sex, across a large number of categories. It is very similar to this example:

```no-highlight
library(brms)
library(palmerpenguins)

penguin <- brm(
  bf(
    body_mass_g ~ 0 + sex + (0 + sex | i | year),
    sigma ~ 0 + sex + (0 + sex | i | year)
  ),
  family = 'lognormal', data = penguins, cores = 4
)

```

(This model doesn’t fit well, but my actual model fits fine.)

Both the mean and the standard deviation depend on sex, and they can vary across the categories (years in this case). In my real data, the mean and `sigma` are strongly related, and this correlation is captured by the `| i |` term.

Now, I want to use this model to get posterior prediction for `sigma`, for a previously unseen category. However, I want this prediction to be conditional on the mean, since this strongly effects the expectation for `sigma`. But, when I run this code, the mean-variance relationship is not taken into account, since the result is the same if I include the mean or not:

```no-highlight
posterior_predict(
  penguin, 
  newdata = data.frame(sex = 'male', year = 2024, body_mass_g = 1e20), 
  allow_new_levels = TRUE,
  dpar = 'sigma'
) |> posterior_summary()

```

The mean-variance correlation is basically zero here, but my actual model it is very clearly not using the response variable, which I think does make sense generally. But this leaves my question:

Is it possible to make use of the mean when predicting sigma?

---

<div class="post-metadata">

**Author:** ![jd\_c](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jd_c/32/16759_2.png) [@jd\_c](https://discourse.mc-stan.org/u/jd_c)\
**Post date:** [October 17, 2024, 6:32pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/2 "2024-10-17T18:32:16Z")

</div>

> [@Ax3man](#):
>
> But, when I run this code, the mean-variance relationship is not taken into account, since the result is the same if I include the mean or not:

`body_mass_g` is the outcome variable… I’m not sure how you are “condition[ing] on the mean” in your predict code here. In your model, I think what you have done is simply correlate `mu[n]` and `sigma[n]` (use the `stancode()` function on your fitted model, or the `make_stancode()` function on the example above) via the varying intercepts and slopes. So, I think when you simply use the `posterior_predict()` function, you are making predictions from parameters that are correlated (I assume. You would need to dig into the details of the `posterior_predict()` function in brms. In any case, I would not think you would get a difference in the results of `posterior_predict()` by including `body_mass_g`, the outcome variable, in the call.)

---

<div class="post-metadata">

**Author:** ![Ax3man](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ax3man/32/11894_2.png) [@Ax3man](https://discourse.mc-stan.org/u/Ax3man)\
**Post date:** [October 17, 2024, 7:02pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/3 "2024-10-17T19:02:10Z")

</div>

Right, I think I understand the model, and I can see why `posterior_predict()` behaves this way.

I am just trying to figure out if I can make a prediction for the `sigma` of a group, and make use of the fact that I have an estimate of `mu` for that group, given that the model includes their correlation among groups. It seems to me that I can make a better prediction for `sigma` if I know `mu`.

---

<div class="post-metadata">

**Author:** ![joels](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/joels/32/10442_2.png) [@joels](https://discourse.mc-stan.org/u/joels)\
**Post date:** [October 17, 2024, 7:07pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/4 "2024-10-17T19:07:27Z")

</div>

Aren’t you already getting the predictive benefit of “knowing `mu`” for each group, in the sense that the posterior fit has already taken into account the correlation between `mu` and `sigma`, since you included that in the model.

Or, to put it another way, if you don’t model `sigma` at all, or model `sigma` without an induced correlation with `mu`, do you get worse predictions than with the model that includes correlations?

---

<div class="post-metadata">

**Author:** ![Ax3man](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ax3man/32/11894_2.png) [@Ax3man](https://discourse.mc-stan.org/u/Ax3man)\
**Post date:** [October 17, 2024, 7:18pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/5 "2024-10-17T19:18:51Z")

</div>

> [@joels](#):
>
> Or, to put it another way, if you don’t model `sigma` at all, or model `sigma` without an induced correlation with `mu`, do you get worse predictions than with the model that includes correlations?

Since I’m interested in `sigma` (not `mu`), I only tried the second. And yes, the model is much better when including their correlation (at least the uncertainty around `sigma` for any particular group is smaller). So in that sense the model is useful.

But I am trying to describe what the expected `sigma` (and the associated uncertainty) is for an unseen group. In my real data the groups are genes, and the outcome is expression. So I want to describe what the expected variance in expression is for a gene. But because this depends strongly on the average gene expression, I would like to describe what the expected expression variance would be, given a certain average expression.

If I just ask for a prediction of `sigma` for an unseen gene, this is very uncertain, since the range of mean expression is very wide.

Hope that makes some sense…

---

<div class="post-metadata">

**Author:** ![joels](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/joels/32/10442_2.png) [@joels](https://discourse.mc-stan.org/u/joels)\
**Post date:** [October 17, 2024, 8:12pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/6 "2024-10-17T20:12:16Z")

</div>

> [@Ax3man](#):
>
> So I want to describe what the expected variance in expression is for a gene. But because this depends strongly on the average gene expression, I would like to describe what the expected expression variance would be, given a certain average expression.

What about using the predicted mean gene expression as a predictor of sigma? (I think this can be done in brms, but might require a nonlinear formula.) I’m not sure if this is kosher, but just thinking out loud.

---

<div class="post-metadata">

**Author:** ![jsocolar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jsocolar/32/2486_2.png) [@jsocolar](https://discourse.mc-stan.org/u/jsocolar)\
**Post date:** [October 18, 2024, 4:53am UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/7 "2024-10-18T04:53:19Z")

</div>

Use `posterior_linpred(... , re_formula = ~ 0, dpar = "sigma")` to get the posterior linpred for `sigma` without the random effect component. Then add in the `sigma` margin of the multivariate random effect conditional on the known `mu` value (note that the distribution for a subset of variables from a multivariate normal, conditional on known values for another subset of variables, is itself a multivariate normal distribution; [this page](https://online.stat.psu.edu/stat505/lesson/6/6.1#:~:text=Conditional%20Distribution%20Properties&text=Any%20distribution%20for%20a%20subset,is%20a%20multivariate%20normal%20distribution.) gives the bivariate case succinctly).

---

<div class="post-metadata">

**Author:** ![joels](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/joels/32/10442_2.png) [@joels](https://discourse.mc-stan.org/u/joels)\
**Post date:** [October 21, 2024, 3:46pm UTC](https://discourse.mc-stan.org/t/predicting-the-variance-taking-the-mean-into-account/36955/8 "2024-10-21T15:46:39Z")

</div>

@Ax3man, can you please post the code you ended up using to implement the approach suggested by jsocolar?