# Model selection with loo\_compare and loo\_model\_weights

**URL:** <https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489>\
**Category:** Modeling\
**Tags:** loo, rstan, techniques\
**Created:** [February 21, 2023, 3:37pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489 "2023-02-21T15:37:17Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [February 21, 2023, 3:37pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/1 "2023-02-21T15:37:17Z")

</div>

I am conducting an analysis of various effects on the harvest of different fish species. Each species is analyzed independent using the same set of 8 nested models. I will not be using the results for predictions, I am trying to capture those predictors which modify the response variable and in what direction.  
It was suggested to me to use both loo\_compare and loo\_model\_weights for model selection. For some of the species, I have come across the situation where the model selected by loo\_compare is not the same as the model with the highest stacking weight by loo\_model\_weights and do not know how best to proceed with reporting results and discussion.

I have read some papers and the various documentation, but it is still unclear to me which method would be best for my study, if another method would be better, or if I should report the results using both methods (although this sounds like it could be confusing for the reader).

---

<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:** [February 21, 2023, 6:25pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/2 "2023-02-21T18:25:26Z")

</div>

Can you show the output you did get from `loo_compare()` and `loo_model_weights()`? I could then provide more detailed answer

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [February 21, 2023, 7:32pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/3 "2023-02-21T19:32:56Z")

</div>

Example 1

```plaintext
loo_compare
      elpd_diff se_diff
mod4 0.0 0.0 
mod10 -5.5 4.4 
mod1 -6.9 6.9 
mod2 -11.8 8.8 
mod3 -26.2 9.8 
mod9 -35.6 13.5 
mod11 -35.9 7.9 
mod5 -399.7 31.3 
mod6 -418.2 31.3 
mod13 -418.3 29.7 
mod8 -418.8 30.2 
mod7 -464.7 34.2 
mod12 -475.8 34.2 
mod14 -572.5 38.2 

loo_model_weights
Method: stacking
------
      weight
mod1 0.180 
mod2 0.106 
mod3 0.202 
mod4 0.097 
mod5 0.000 
mod6 0.000 
mod7 0.001 
mod8 0.000 
mod9 0.000 
mod10 0.191 
mod11 0.091 
mod12 0.000 
mod13 0.106 
mod14 0.025 

```

Example 2

```plaintext
loo_compare
> loo_out
      elpd_diff se_diff
mod3 0.0 0.0 
mod1 -0.4 3.4 
mod11 -1.7 6.3 
mod4 -2.8 5.7 
mod10 -91.1 14.0 
mod2 -94.2 13.6 
mod9 -96.3 13.5 
mod5 -153.9 16.1 
mod8 -175.5 19.1 
mod6 -253.7 19.1 
mod13 -278.0 21.0 
mod7 -346.4 28.2 
mod14 -398.9 24.4 
mod12 -402.2 29.5 

 loo_model_weights
Method: stacking
------
      weight
mod1 0.148 
mod2 0.000 
mod3 0.000 
mod4 0.011 
mod5 0.064 
mod6 0.000 
mod7 0.137 
mod8 0.122 
mod9 0.000 
mod10 0.194 
mod11 0.323 
mod12 0.002 
mod13 0.000 
mod14 0.000 

```

I should add… that these are not exactly nested models because I checked each model for multicollinearity and removed the highest predictor with a VIF\>5, this process was then repeated to bring all VIF \< 5.

---

<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:** [February 21, 2023, 8:42pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/4 "2023-02-21T20:42:19Z")

</div>

> [@Jacob](#):
>
> ```no-highlight
> loo_compare
> elpd_diff se_diff
> mod4 0.0 0.0 
> mod10 -5.5 4.4 
> mod1 -6.9 6.9 
> mod2 -11.8 8.8 
> 
> ```

These are so close to each other that we can’t say that one would be better than others. See more information at [CV-FAQ 16](https://avehtari.github.io/modelselection/CV-FAQ.html#se_diff) and the reference given there.

Stacking is not model selection, but model averaging and can give non-zero weight for a model that is on average much worse, but the mixture of weighted models can still be better. For example if the true distribution would be t-distribution and you have one narrow and one wide normal, one of the normal’s can be much worse alone, but the mixture can be closer to t-distribution than either normal, and thus the weights can be non-zero for both.

Your stacking weights are mostly non-zero for the models that have good predictive performance, but for example mod13 has weight 0.106, although it’s predictive performance alone is much worse than the performance of the best ones. This is likely due to mod13 having thicker tails in the predictive distribution and adding that to the mixture improves the combination just a little. This hints that it might be possible to modify one of the top models to have an observation model that better matches the true distribution.

> [@Jacob](#):
>
> Example 2
> 
> ```no-highlight
> loo_compare
> > loo_out
> elpd_diff se_diff
> mod3 0.0 0.0 
> mod1 -0.4 3.4 
> mod11 -1.7 6.3 
> mod4 -2.8 5.7 
> 
> ```

These models are indistinguishable based on the predictive performance and the rest have much worse predictive performance.

> [@Jacob](#):
>
> I should add… that these are not exactly nested models because

I didn’t assume they are nested, and the stacking weights even hint that they are not nested as so many models have non-zero weight.

> [@Jacob](#):
>
> because I checked each model for multicollinearity and removed the highest predictor with a VIF\>5, this process was then repeated to bring all VIF \< 5.

This will make the model comparison very difficult. Instead, I recommend using all the covariates and use projection predictive variable selection, especially if the only differences between the models are which covariates are used. If you tell more about the models, I may be able to provide further advice.

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [February 21, 2023, 9:30pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/5 "2023-02-21T21:30:44Z")

</div>

Thank you so much for engaging with me. What else would you like to know about the models?  
I will read your article on Projection predictive variable selection.

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [February 22, 2023, 3:34pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/6 "2023-02-22T15:34:51Z")

</div>

I am fitting to a neg\_binomial\_2 family model, which does not seem to be supported by projpred. Here are two figures, first the ppc\_dens\_overlay, and second the mcmc\_areas for the full model for one of my species (one with a particularly large number of 0 estimates!).

 ![ppc_dens_overlay](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/3e10f14240d0c15a41b2e3ae4035476b048e78e0.png)  
 ![mcmc_areas](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/2/231313d7dd69b927205b2dc042993308ceebaad0.png)

---

<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:** [February 27, 2023, 7:35pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/7 "2023-02-27T19:35:03Z")

</div>

> [@Jacob](#):
>
> I am fitting to a neg\_binomial\_2 family model, which does not seem to be supported by projpred.

I think you could do this using the recently added latent space approach. The vignette has also negative binomial example [Latent projection predictive feature selection](https://cran.r-project.org/web/packages/projpred/vignettes/latent.html)

The marginal posterior hints collinearity, and projpred would probably work well

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [March 3, 2023, 3:34pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/8 "2023-03-03T15:34:39Z")

</div>

Thank you.  
I got as far as identifying the solution terms ( `soltrms_lat_final` in the example).  
The example then refers to the main vignette for post-selection inference. That vignette calls the function `project` which does not accept neg\_binomial\_2 models.

I do not need to make predictions from the final model, but am rather looking for influential variables and the direction they modify the response variable. Can I simply use the solution terms identified to rerun the model using the whole data set?

---

<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 4, 2023, 4:49pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/9 "2023-03-04T16:49:44Z")

</div>

Pinging @fweber144 who made the vignette, maybe he is able to help

---

<div class="post-metadata">

**Author:** ![fweber144](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fweber144/32/15592_2.png) [@fweber144](https://discourse.mc-stan.org/u/fweber144)\
**Post date:** [March 4, 2023, 7:48pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/10 "2023-03-04T19:48:20Z")

</div>

> [@Jacob](#):
>
> The example then refers to the main vignette for post-selection inference. That vignette calls the function `project` which does not accept neg\_binomial\_2 models.

Did you set `latent = TRUE` and the other arguments relevant for the latent projection (here: `latent_ll_oscale = latent_ll_oscale_nebin` and `latent_ppd_oscale = latent_ppd_oscale_nebin`) also in the `project()` call, not only in the `varsel()` or `cv_varsel()` call? This is explained in section [Implementation](https://mc-stan.org/projpred/articles/latent.html#implementation), but I guess I should have pointed that out more clearly.

> [@Jacob](#):
>
> I do not need to make predictions from the final model, but am rather looking for influential variables and the direction they modify the response variable. Can I simply use the solution terms identified to rerun the model using the whole data set?

That sounds like you don’t have an actual prediction problem but rather want to interpret the parameter estimates. Just be warned that projpred is not made for parameter interpretation, but prediction problems.

---

<div class="post-metadata">

**Author:** ![fweber144](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fweber144/32/15592_2.png) [@fweber144](https://discourse.mc-stan.org/u/fweber144)\
**Post date:** [March 4, 2023, 8:37pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/11 "2023-03-04T20:37:35Z")

</div>

You can avoid the repetition of the arguments related to the latent projection by defining the reference model object of class `refmodel` first (via `projpred::get_refmodel()` which redirects to `projpred:::get_refmodel.stanreg()` in your case; it is in this creation of the reference model object where you add the arguments related to the latent projection) and then calling `varsel()`, `cv_varsel()`, or `project()` on that predefined reference model object (instead of on the reference model _fit_) without the need for repeating the arguments related to the latent projection.

Alternatively, if you have run the code from the latent-projection vignette, you can call `project()` on the output of `varsel()` or `cv_varsel()` (instead of on the reference model _fit_) without the need for repeating the arguments related to the latent projection.

All of this might sound a bit complicated now but is in fact a quite natural approach when you get more familiar with the UI of projpred. Anyway, I guess I really need to revise that part of the latent-projection vignette. Sorry for that. I’ll put it on the issue tracker.

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [March 9, 2023, 11:27am UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/13 "2023-03-09T11:27:35Z")

</div>

Again, thank you both for your help.  
I was able to get the projections to work using the latent = TRUE argument.  
Getting back to my problem. @avehtar I started down this path at your suggestion that I can use this approach instead of model selection. But once I have selected the best predictor terms in my model, I need to deviate from the vignettes to use those terms for parameter interpretation. Is it simply enough to rerun the model using the terms identified?

---

<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 20, 2023, 1:47pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/14 "2023-03-20T13:47:21Z")

</div>

> [@Jacob](#):
>
> Getting back to my problem. @avehtar I started down this path at your suggestion that I can use this approach instead of model selection. But once I have selected the best predictor terms in my model, I need to deviate from the vignettes to use those terms for parameter interpretation. Is it simply enough to rerun the model using the terms identified?

I missed this earlier (you tried to tag me, but missed one letter in my username).

1. For the parameter interpretation you need to take into account causal assumptions and possible collinearities
2. If you want to be careful, we recommend doing the interpretation with the reference model (the big model). Although the projection has nice decision theoretical interpretation, the actual implementation is such that there can be difference to the theory and we don’t have yet such diagnostic that we would be happy recommending using the projected posterior for interpretation (assuming causal assumptions allow that). For posterior interpretation (assuming causal assumptions allow that) the refitted model posterior with the selected variables, is likely to be safer than the projected posterior (the variable selection procedure itself is stable), but still we recommend to make possible causal interpretations using the reference model.

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [April 3, 2023, 1:37pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/15 "2023-04-03T13:37:16Z")

</div>

@ [avehtari](https://discourse.mc-stan.org/u/avehtari)

Thanks for the tips. This approach seems to work very well for most of the species I am analyzing, but a few species have issues with the reference model. This is particularly the case for the species with fewer observations.

For example, for this species the mlpd plot:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/3/8/380bb1ebb450995123c250acbee9bc0f234ce090.png)  
and with ‘deltas = FALSE’ :  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/d/5/d58bdfa1c18f045b54048a336601ff84fc3772a5.png)  
suggesting that the reference model is not good. Are there other ways I can improve the reference model? (I can’t add more observations). Here is the code for the reference model (following the Latent projection predictive feature selection example)

D ← 12  
p0 ← 6  
mu\_prior ← 100  
tau0 ← p0 / (D - p0) / sqrt(mu\_prior) / sqrt(N)

refm\_fit\_nebin ← stan\_glm(  
formula = har ~  
p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + p10 + p11 + p12,  
family = neg\_binomial\_2(),  
data = train,  
prior = hs(global\_scale = tau0, slab\_df = 100, slab\_scale = 1),  
seed = 7304, QR = TRUE,  
adapt\_delta = 0.99,  
iter = 15000,  
cores = 4  
)

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [April 3, 2023, 3:44pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/16 "2023-04-03T15:44:30Z")

</div>

@fweber144

I just saw this note “For the sake of simplicity, we won’t adjust `tau0` to this new family, but in a real-world example, such an adjustment would be necessary. However, since Table 1 of Piironen and Vehtari ([2017](https://cran.r-project.org/web/packages/projpred/vignettes/latent.html#ref-piironen_sparsity_2017)) doesn’t list the negative binomial distribution, this would first require a manual derivation of the pseudo-variance σ~2”

Has pseudo-variance for negative binomial distributions been estimated since you wrote this article?

---

<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, 2023, 5:08pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/17 "2023-04-03T17:08:28Z")

</div>

> [@Jacob](#):
>
> suggesting that the reference model is not good

I interpret the plots so that reference model is useful. Null model is much worse, and you need many variables before reaching the reference model performance.

---

<div class="post-metadata">

**Author:** ![fweber144](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fweber144/32/15592_2.png) [@fweber144](https://discourse.mc-stan.org/u/fweber144)\
**Post date:** [April 3, 2023, 7:59pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/18 "2023-04-03T19:59:17Z")

</div>

> Has pseudo-variance for negative binomial distributions been estimated since you wrote this article?

No, unfortunately not, at least to my knowledge. If you don’t want to derive it on your own (or if it’s not possible; I haven’t tried), you could try to use the pseudo-variance for the Poisson distribution (third line in Table 1 of Piironen and Vehtari, [2017](https://cran.r-project.org/web/packages/projpred/vignettes/latent.html#ref-piironen_sparsity_2017)). But that would be just an approximation.

---

<div class="post-metadata">

**Author:** ![Jacob](https://avatars.discourse-cdn.com/v4/letter/j/f1d935/32.png) [@Jacob](https://discourse.mc-stan.org/u/Jacob)\
**Post date:** [July 14, 2023, 1:43pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/19 "2023-07-14T13:43:07Z")

</div>

Looping back to this comment, I have some clarifying questions.

> [@avehtari](#):
>
> If you want to be careful, we recommend doing the interpretation with the reference model (the big model).

I assume by ‘the reference model (the big model)’ you mean the model with all of the predictors (for which the vignette uses the training data). To clarify, I randomly sampled from my observations so that half of the observations were used as the training data and half as the test data.

> [@avehtari](#):
>
> For posterior interpretation (assuming causal assumptions allow that) the refitted model posterior with the selected variables, is likely to be safer than the projected posterior (the variable selection procedure itself is stable), but still we recommend to make possible causal interpretations using the reference model.

Are you saying that for posterior interpretation refitting the model with the selected variables is safer than projecting the posterior? If so, do you use all of the observations for refitting with the selected variables or do you use the test data? I am having a hard time here understanding how the refitted model posterior with the selected variables can be used to make causal interpretations using the reference model if the reference model itself has all of the variables rather than the selected variables.

---

<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:** [July 25, 2023, 12:10pm UTC](https://discourse.mc-stan.org/t/model-selection-with-loo-compare-and-loo-model-weights/30489/20 "2023-07-25T12:10:41Z")

</div>

> [@Jacob](#):
>
> Are you saying that for posterior interpretation refitting the model with the selected variables is safer than projecting the posterior?

Yes. Then we need to worry only about the selection induced bias (which is reduced by projection predictive selection), but no need to worry the error from the approximate projection. There is some discussion in [[2306.15581] Robust and efficient projection predictive inference](https://arxiv.org/abs/2306.15581)

> [@Jacob](#):
>
> If so, do you use all of the observations for refitting with the selected variables or do you use the test data?

I always use all the data with combination of using cross-validation to assess the out-of-sample predictive performance. See also [[2306.15581] Robust and efficient projection predictive inference](https://arxiv.org/abs/2306.15581)

> [@Jacob](#):
>
> I am having a hard time here understanding how the refitted model posterior with the selected variables can be used to make causal interpretations using the reference model if the reference model itself has all of the variables rather than the selected variables.

We primarily recommend doing the causal interpretation using the big model which has been constructed taking into account the causal constraints. Variable selection can be done to find irrelevant variables, but depending on the causal structure, the variable selection is valid only for some covariates. If you have made the variable selection honoring the causal constraints, then you could use the smaller model for the causal interpretation, but there is no benefit from using the smaller model for the causal interpretation.