# Calculating LOO-CV for a multinormal regression model

**URL:** <https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483>\
**Category:** Modeling\
**Tags:** loo\
**Created:** [March 3, 2020, 10:06am UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483 "2020-03-03T10:06:46Z")\
**Posts on this page:** 13\
**Page:** 2

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 25, 2020, 8:19pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/21 "2020-03-25T20:19:25Z")

</div>

> [@ofirshan](#):
>
> > Can you compute also pointwise lpd?
> 
> I guess this means to run loo on the log\_lik matrix I extracted from the k-fold?

No. lpd is defined in equation (2) in [[1507.04544] Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC](https://arxiv.org/abs/1507.04544)

```no-highlight
lpd <- matrixStats::colLogSumExps(ll) - log(nrow(ll))

```

where `ll` is log\_lik from the full data posterior. Then you can compute pointwise p\_loo by subtracting elpd\_kcv from lpd. If you then plot observation index vs pointwise p\_loo, we’ll learn more. I’m sorry my answers are brief as I’m heavily overloaded, but I still want to help.

---

<div class="post-metadata">

**Author:** ![ofirshan](https://avatars.discourse-cdn.com/v4/letter/o/cc9497/32.png) [@ofirshan](https://discourse.mc-stan.org/u/ofirshan)\
**Post date:** [March 26, 2020, 1:26pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/22 "2020-03-26T13:26:33Z")

</div>

Hi,

So I hope I got it right: I ran the equation on both the log-lik from the full data posterior (this is lpd), as well as on the log-lik from the heldout data posterior (this is elpd\_kcv). Substracting these two yields an _n_-long, 2000 pt. vector, which is the pointwise p\_loo - correct? As you suggested, I plotted it against the observations indices (I’m attaching two types of visualizations. It’s still very dense, but hopefully it provides you with sufficient information; if not’ please let me know what would be a more suitable visualization:

 ![pw_pLoo](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/12543fa0129040ef45e038249f2ca37b5ec294e1.png)

Thanks again for all you tremendous help,

Ofir

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 26, 2020, 4:12pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/23 "2020-03-26T16:12:58Z")

</div>

> [@ofirshan](#):
>
> which is the pointwise p\_loo - correct?

Yes. And the plot looks sensible. For models with conditionally univariate observation model pointwise p\_loo is usually \< 1 unless there is really bad model misspecification. With multinormal observation model this can be more, but I haven’t had experience with than and haven’t had time to go through the math to figure out the asymptotic behavior in weak prior well specified model case.

So it seems that when you remove one observation the posterior is changing too much for PSIS to work. But since you got KCV to work you are fine. The model is very flexible. You can use KCV also with other utility or loss function. What do you think would be useful measure for accuracy in your case?

---

<div class="post-metadata">

**Author:** ![ofirshan](https://avatars.discourse-cdn.com/v4/letter/o/cc9497/32.png) [@ofirshan](https://discourse.mc-stan.org/u/ofirshan)\
**Post date:** [March 28, 2020, 7:00pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/24 "2020-03-28T19:00:08Z")

</div>

> What do you think would be useful measure for accuracy in your case?

I thought that this is somehow indicated by the se\_elpd, and that the next and final step would be to compare this between models (at least that was the procedure described in the tutorial I followed).  
Do I actually need to generate data within each fold by using the holdout data and the estimated model parameters, via the **`multi_normal_cholesky_rng`** function? and then maybe compare the standard deviation of this simulated data vs. that of the true data? And if this is the case, should I do it for each model before or after model comparison (is it possible for one model to outperform another based on se\_elpd comparison, yet to be less accurate?)? Is there as simpler solution or common practice i’m missing here?

Thanks,

Ofir

---

<div class="post-metadata">

**Author:** ![ofirshan](https://avatars.discourse-cdn.com/v4/letter/o/cc9497/32.png) [@ofirshan](https://discourse.mc-stan.org/u/ofirshan)\
**Post date:** [March 28, 2020, 9:13pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/25 "2020-03-28T21:13:31Z")

</div>

> and then maybe compare the standard deviation of this simulated data vs. that of the true data?

- this is a mistake, I meant here something like averaging the squared error between the predicted heldout data and the true heldout data.

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 31, 2020, 6:28pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/26 "2020-03-31T18:28:00Z")

</div>

> [@ofirshan](#):
>
> I thought that this is somehow indicated by the se\_elpd

se\_elpd estimates the accuracy of elpd\_loo, but what I meant by my question is that what do you think is useful measure for accuracy of your model? That, is if you make a prediction with your model how would you like to compare that to the actual observation? What kind of comparison would be useful, meaningful and interpretable taking into account your modeling problem?

> [@ofirshan](#):
>
> indicated by the se\_elpd, and that the next and final step would be to compare this between models (at least that was the procedure described in the tutorial I followed).

Maybe you have type there as you would not compare SE’s which measure the accuracy of the estimate.

> [@ofirshan](#):
>
> Do I actually need to generate data within each fold by using the holdout data and the estimated model parameters, via the **`multi_normal_cholesky_rng`** function?

It depends what utility or loss function you want to use.

> [@ofirshan](#):
>
> and then maybe compare the standard deviation of this simulated data vs. that of the true data?

That goes in to the direction of posterior predictive checking which is useful, too. See e.g. [Visualization in Bayesian workflow](https://doi.org/10.1111/rssa.12378)

> [@ofirshan](#):
>
> And if this is the case, should I do it for each model before or after model comparison (is it possible for one model to outperform another based on se\_elpd comparison, yet to be less accurate?)?

Again that se\_elpd in wrong context?

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 31, 2020, 6:29pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/27 "2020-03-31T18:29:10Z")

</div>

> [@ofirshan](#):
>
> - this is a mistake, I meant here something like averaging the squared error between the predicted heldout data and the true heldout data.

Sure, if you think that is easy for you to interpret and it’s meaningful in your modeling problem.

---

<div class="post-metadata">

**Author:** ![ofirshan](https://avatars.discourse-cdn.com/v4/letter/o/cc9497/32.png) [@ofirshan](https://discourse.mc-stan.org/u/ofirshan)\
**Post date:** [March 31, 2020, 6:45pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/28 "2020-03-31T18:45:31Z")

</div>

Hi,  
I indeed had a typo their with the se\_elpd.  
I finished comparing the 3 models. I did this by creating a list for each model, which contains the pointwise elpd estimate and se. I constructed the structures so I could run loo\_comapre on these lists. Running loo\_compare(model\_1, model\_2, model\_3) worked fine - is this the right way to compare? Because if I added (criterion=“kfold”) I got an error saying “All inputs should have class ‘loo’”; I don’t understand this, since I did define class ‘loo’ and ‘kfold’.  
So these are my results:

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/02f87fad068aa959e9c0d37f32706d50db0f2225.png)

Are these differences considered “significant”? I thought that it’s a bit strange that model 2 has a much smaller se. Anyway, i’m now running it again with different holdout indices, just to be sure.

Many thanks again,

Ofir

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [April 1, 2020, 7:37pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/29 "2020-04-01T19:37:08Z")

</div>

> [@ofirshan](#):
>
> I indeed had a typo their with the se\_elpd.

What did you mean to write?

> [@ofirshan](#):
>
> I finished comparing the 3 models. I did this by creating a list for each model, which contains the pointwise elpd estimate and se.

How did you estimate these? The thread is long and you have reported different results, so I’m uncertain what computations you have used for these. You can post the code lines you are using so it’s easier to check what you’ve done.

> [@ofirshan](#):
>
> Running loo\_compare(model\_1, model\_2, model\_3)

And what are these model\_x objects?

> [@ofirshan](#):
>
> Because if I added (criterion=“kfold”)

Added where?

> [@ofirshan](#):
>
> Are these differences considered “significant”?

Are these based on K-fold-CV? And is the table from loo\_compare?

> [@ofirshan](#):
>
> I thought that it’s a bit strange that model 2 has a much smaller se.

It’s not much smaller and different models may have different variation in predictive accuracy affecting the accuracy of expectation.

---

<div class="post-metadata">

**Author:** ![ofirshan](https://avatars.discourse-cdn.com/v4/letter/o/cc9497/32.png) [@ofirshan](https://discourse.mc-stan.org/u/ofirshan)\
**Post date:** [April 1, 2020, 8:43pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/30 "2020-04-01T20:43:41Z")

</div>

Hi,

> What did you mean to write?

I meant to write elpd estimate.

> How did you estimate these? The thread is long and you have reported different results, so I’m uncertain what computations you have used for these. You can post the code lines you are using so it’s easier to check what you’ve done.

I ran k-fold with 10 folds for each of my three models in rstan. Then I extracted log-lik and computed pointwise elpd for the k-fold, based on the functions in this link - [K-fold cross-validation in Stan | DataScience+](https://datascienceplus.com/k-fold-cross-validation-in-stan/). So this was the function for extracting elpd estimate and its se:

> kfold ← function(log\_lik\_heldout) {  
> library(matrixStats)  
> logColMeansExp ← function(x) {  
> # should be more stable than log(colMeans(exp(x)))  
> S ← nrow(x)  
> colLogSumExps(x) - log(S)  
> }

# See equation (20) of @VehtariEtAl2016

pointwise ← matrix(logColMeansExp(log\_lik\_heldout), ncol= 1)  
colnames(pointwise) ← “elpd”

# See equation (21) of @VehtariEtAl2016

elpd\_kfold ← sum(pointwise)  
se\_elpd\_kfold ← sqrt(ncol(log\_lik\_heldout) \* var(pointwise))  
out ← list(  
pointwise = pointwise,  
elpd\_kfold = elpd\_kfold,  
se\_elpd\_kfold = se\_elpd\_kfold)  
#structure(out, class = “loo”)  
return(out)  
}

> And what are these model\_x objects?

Based on what you wrote in this thread - [Compare models with k-fold cv - #2 by avehtari](https://discourse.mc-stan.org/t/compare-models-with-k-fold-cv/9042/2), for each model I created a structure consisting of the pointwise elpd estimate and its se derived from the kfold. The structure was defined according to the guidelines in the loo package manual regarding kfold (pp. 12, under “kfold generic”, so it will be compatible with other loo functions - hope I got it right.

> Added where?

In the loo\_comapre function, I thought I’m supposed to define this if i’m performing it on k-fold.

> Are these based on K-fold-CV? And is the table from loo\_compare?

Yes, based on th kcv, and this is the table yielded by the loo-compare.

Thanks,

Ofir

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [April 3, 2020, 6:44pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/31 "2020-04-03T18:44:05Z")

</div>

> [@ofirshan](#):
>
> In the loo\_comapre function, I thought I’m supposed to define this if i’m performing it on k-fold.

loo\_compare doesn’t have that kind of argument [Model comparison — loo\_compare • loo](https://mc-stan.org/loo/reference/loo_compare.html), so I’m confused.

> [@ofirshan](#):
>
> Yes, based on th kcv, and this is the table yielded by the loo-compare.

Good. The values seem reasonable given the information you have provided.

---

<div class="post-metadata">

**Author:** ![ofirshan](https://avatars.discourse-cdn.com/v4/letter/o/cc9497/32.png) [@ofirshan](https://discourse.mc-stan.org/u/ofirshan)\
**Post date:** [April 6, 2020, 1:19pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/32 "2020-04-06T13:19:21Z")

</div>

> loo\_compare doesn’t have that kind of argument [https://mc-stan.org/loo/reference/loo\_compare.html](https://mc-stan.org/loo/reference/loo_compare.html), so I’m confused.

I saw it here - [Information criteria and cross-validation — loo.stanreg • rstanarm](http://mc-stan.org/rstanarm/reference/loo.stanreg.html) - maybe it’s only relevant for fits derived from rstanarm?

> Good. The values seem reasonable given the information you have provided.

Great! and thank you so much for all your help along this process.

Best regards,

Ofir

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [April 7, 2020, 6:56pm UTC](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483/33 "2020-04-07T18:56:55Z")

</div>

> [@ofirshan](#):
>
> I saw it here - [http://mc-stan.org/rstanarm/reference/loo.stanreg.html](http://mc-stan.org/rstanarm/reference/loo.stanreg.html) - maybe it’s only relevant for fits derived from rstanarm?

It’s relevant when the used for stanreg objects from rstanarm for which we can make automated kfold, but we can’t make automated kfold for stanfit objects. I’ll make a not to clarify that function reference.

[Previous page](https://discourse.mc-stan.org/t/calculating-loo-cv-for-a-multinormal-regression-model/13483.md?page=1)