# Projpred elpd goes down and RMSE goes up after x variables

**URL:** <https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153>\
**Category:** Modeling\
**Tags:** loo\
**Created:** [February 18, 2020, 3:55am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153 "2020-02-18T03:55:27Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![abartonicek](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abartonicek/32/7020_2.png) [@abartonicek](https://discourse.mc-stan.org/u/abartonicek)\
**Post date:** [February 18, 2020, 3:55am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/1 "2020-02-18T03:55:27Z")

</div>

Hi, I’m using the `projpred` to find a smaller submodel that predicts similarly to a full regularized horseshoe model with 31 predictors. I used the `cv_varsel()` function to fit the submodels & compare them using PSIS-LOO. Now, when I use the `suggest_size()` function, it suggests 12 variables as the optimal number of variables to approximate the predictions of the full model. However, when I plot the trajectory of the selection using the `varsel_plot()` function, I get the following plot:

 ![projpred](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/4/4e9bd1225e44d8c38faeef6c6ed363aa552fd9f0.png)

Does the model start to overfit after 6 predictors & give worse out-of-sample predictions as a result (and therefore the optimal number of variables in the submodel would actually be 6)? Or is there something wrong-er going on?

---

<div class="post-metadata">

**Author:** ![mcol](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/mcol/32/5598_2.png) [@mcol](https://discourse.mc-stan.org/u/mcol)\
**Post date:** [February 18, 2020, 7:05am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/2 "2020-02-18T07:05:02Z")

</div>

Yes, I would interpret this as overfitting. I’ve never found `suggest_size` to be particularly helpful, as it returns only a number, and often, as in your case, it can pick a more complex model than perhaps necessary. The graph is much clearer and lets you judge better the trade-off betweem model size and predictive performance.

---

<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 18, 2020, 8:02am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/3 "2020-02-18T08:02:55Z")

</div>

I would not interpret this as overfitting. Overfitting could be part of the problem, but that term oversimplifies the issue and it’s not even well defined term (or if you have a good definition please tell).

My guess is correlating covariates with outliers in covariate values, and the model is misspecified and not modeling them appropriately. `varsel_plot` shows cross-validation results, and cross-validation can have some bias in case of misspecified models. I recommend to check correlation plot for the covariates and in more detail check what covariates are in the group from 7-12 and if they have suspicious values.

We use `suggest_size()` mostly when we need to run 1000 simulation experiments and want to automatically select the size. Otherwise we highly recommend to look at the plots as there might be implicit preference for a smaller model and the user can then decide how much predictive performance they are willing to sacrifice. The current `varsel_plot()` and `suggest_size()` are suffering a bit from the behavior of cross-validation when comparing models which are very similar. We are in progress of writing a paper describing the issue in more detail, and we have a better idea for `suggest_size()` which will hopefully be more robust.

I may be able to provide additional comments if you tell also

- `loo()` output for the dull model
- did you normalize covariates and target?
- how did you select the regularized horseshoe prior parameter?

---

<div class="post-metadata">

**Author:** ![abartonicek](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abartonicek/32/7020_2.png) [@abartonicek](https://discourse.mc-stan.org/u/abartonicek)\
**Post date:** [February 18, 2020, 8:58pm UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/4 "2020-02-18T20:58:47Z")

</div>

Thank you for responding so quickly Aki!

So, if I understand correctly, you suggest that the weird behavior may be the result of collinearity? I’ve checked bivariate correlations between the problematic predictors & there’s few moderately high ones, with the highest one being 0.63. Could that be it? If so, is there a way to address that, perhaps with more regularization? In terms of suspicious values, I’ve double checked the model summary and `mcmc_areas()` plot & I can’t see anything unreasonable, as in there aren’t any slopes much larger than I would expect based on background knowledge (is that what you meant)?

As for your questions:

- I did scale and center all continuous predictors
- I set 1 degree of freedom for the local student-t shrinkage of the horseshoe prior, for everything else I left the default
- Here’s the `loo()` output for the full model:

```
Computed from 8000 by 593 log-likelihood matrix

         Estimate SE
elpd_loo -626.6 19.6
p_loo 24.5 2.0
looic 1253.1 39.1
------
Monte Carlo SE of elpd_loo is 0.1.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.

```

---

<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 19, 2020, 9:15am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/5 "2020-02-19T09:15:14Z")

</div>

> [@abartonicek](#):
>
> So, if I understand correctly, you suggest that the weird behavior may be the result of collinearity?

Not just the collinearity. Look also how the marginal distribution of the covariate values look like. If there are strange values, then including one of the correlating variables can throw the predictions off, but when you include many correlating values then effect of the strange values in each individual covariate will go down.

> [@abartonicek](#):
>
> If so, is there a way to address that, perhaps with more regularization?

Not more regularization, but more modeling. If it’s about strange measurement errors in covariates, you should include measurement error model. There can be also cases, where you just need to include group of correlating covariates for better performance.

> [@abartonicek](#):
>
> In terms of suspicious values, I’ve double checked the model summary

No, look at the covariate values. Are there values which are marginally or jointly far away from the other values?

> [@abartonicek](#):
>
> I set 1 degree of freedom for the local student-t shrinkage of the horseshoe prior, for everything else I left the default

You should set the global scale parameter as instructed in [Sparsity information and regularization in the horseshoe and other shrinkage priors](https://projecteuclid.org/euclid.ejs/1513306866)  
For example, if you assume a priori that there could 5 practically non-zero co-efficients then use

```no-highlight
p0 <- 5 # prior guess for the number of relevant variables
tau0 <- p0/(p-p0) * 1/sqrt(n)
rhs_prior <- hs(global_scale=tau0)

```

where `p` is the total number of covariates and `n` is the number of observations. With your `p` and `n` and `p0=5`, the default global scale is not far away from what the above formula gives.

> [@abartonicek](#):
>
> Here’s the `loo()` output for the full model:

`p_loo` value combined with the information that there are correlating covariates, indicates that most covariates are relevant, so 12 covariates from projpred is not that surprising.

As all Pareto k estimates \< 0.5, there are no really bad outliers, but then you have 20 times more observations than covariates which makes the posterior easier. projpred might work equally well with Gaussian prior.

---

<div class="post-metadata">

**Author:** ![abartonicek](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abartonicek/32/7020_2.png) [@abartonicek](https://discourse.mc-stan.org/u/abartonicek)\
**Post date:** [February 19, 2020, 11:53pm UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/6 "2020-02-19T23:53:23Z")

</div>

Thank you, so I’ve done couple things:

> [@avehtari](#):
>
> Look also how the marginal distribution of the covariate values look like. If there are strange values, then including one of the correlating variables can throw the predictions off, but when you include many correlating values then effect of the strange values in each individual covariate will go down.

I’m sorry I’m still not very well-versed in the terminology, but I’ve looked at the marginal posterior distributions of all predictors individually, and at the pairwise joint marginal distributions as well. The marginal distributions for single predictors look reasonable to me, and I can’t find any joint pairwise posteriors that would not overlap zero (while individual predictors do). Is that the right way to go about it?

> [@avehtari](#):
>
> As all Pareto k estimates \< 0.5, there are no really bad outliers, but then you have 20 times more observations than covariates which makes the posterior easier. projpred might work equally well with Gaussian prior.

I’ve tried to fit the model with a normal(0, 1) prior (a “Bayesian ridge”, the way I understand it), and it seems to give a little more reasonable looking selection trajectory plot:

 ![projpred2](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/c/c744aabd0c28baef776b8743c3a1c0238b2a896d.png)

Additionally, the `suggest_size()` output is now 5. However, `cv_varsel()` now gives the following message:

```no-highlight
Warning messages:
1: posterior_linpred(transform = TRUE) is deprecated.Please use pp_expect() instead. 
2: Quick-TRANSfer stage steps exceeded maximum (= 400000) 
3: Quick-TRANSfer stage steps exceeded maximum (= 400000) 

```

Is that something to worry about?

Also, I have a test set that I can test the submodel predictions on, is that something I should do now or should I wait before I get a more well-behaved feature selection?

Thank you so much for your help

---

<div class="post-metadata">

**Author:** ![abartonicek](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abartonicek/32/7020_2.png) [@abartonicek](https://discourse.mc-stan.org/u/abartonicek)\
**Post date:** [February 20, 2020, 1:38am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/7 "2020-02-20T01:38:04Z")

</div>

Also, is there a way to get `bayes_R2()` for the projected submodel?

EDIT: Computed it using the script in: [https://github.com/jgabry/bayes\_R2/blob/master/bayes\_R2.pdf](https://github.com/jgabry/bayes_R2/blob/master/bayes_R2.pdf)

---

<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 20, 2020, 9:45am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/8 "2020-02-20T09:45:31Z")

</div>

> [@abartonicek](#):
>
> I’m sorry I’m still not very well-versed in the terminology,

No problem, I’ll keep trying to use more explicit sentences.

> [@abartonicek](#):
>
> but I’ve looked at the marginal posterior distributions of all predictors individually, and at the pairwise joint marginal distributions as well.

Don’t look at the posterior. I asked you to look at the data. The data is the one you have before using Stan. I said “covariate values”, you used in the topic title “variable” instead of “covariate”, but that can lead to confusion as the model can also have “variables”. In the above you use “predictor” which is same “covariate” and good word for that, too. The predictor is data and that doesn’t have posterior. The model has coefficients or parameters (or weights) corresponding to the predictors. These coefficients/parameters have posterior. Now look at the data and not the posterior.

> [@abartonicek](#):
>
> it seems to give a little more reasonable looking selection trajectory plot:

This can happen since n\>\>p and the predictors are correlating.

It’s difficult to see from the plot whether the full model performance is similar with RHS or normal prior?

> [@abartonicek](#):
>
> Is that something to worry about?

First warning not. 2. 3. warning I don’t know, but @AlejandroCatalina might know the reason.

> [@abartonicek](#):
>
> Also, I have a test set that I can test the submodel predictions on, is that something I should do now or should I wait before I get a more well-behaved feature selection?

I think they are reasonably well behaved. In both cases there is similar bump in the curve, and it could bee good to check those predictor values.

---

<div class="post-metadata">

**Author:** ![AlejandroCatalina](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/alejandrocatalina/32/6093_2.png) [@AlejandroCatalina](https://discourse.mc-stan.org/u/AlejandroCatalina)\
**Post date:** [February 20, 2020, 9:54am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/9 "2020-02-20T09:54:16Z")

</div>

Regarding the warnings of the form:

> [@abartonicek](#):
>
> Quick-TRANSfer stage steps exceeded maximum (= 400000)

These happen when clustering the data points with `kmeans` and some of the data points are so close together that the algorithm does not converge in the maximum limit of steps. I don’t think these are a reason to worry in this particular case.

Nonetheless, I’ll keep that in mind to make the `kmeans` step a bit more robust now that we are working on improving `projpred` in general!

---

<div class="post-metadata">

**Author:** ![abartonicek](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abartonicek/32/7020_2.png) [@abartonicek](https://discourse.mc-stan.org/u/abartonicek)\
**Post date:** [February 20, 2020, 10:13pm UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/10 "2020-02-20T22:13:07Z")

</div>

> [@avehtari](#):
>
> Now look at the data and not the posterior.

I’ve double checked the data. I noticed there was one continuous covariate that I missed in my data pre-processing pipeline and didn’t standardize (very stupid mistake). I’ve fixed it and now the selection trajectory looks like this:

 ![projpred3](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/b19bc48709dea44501e7c42577227d67cc8ece53.png)

The trajectory looks smoother, there’s still the same bump but it’s smaller now. The `suggest_size()` output is still 5. The unstandardized covariate did not appear until late in the feature selection trajectory, and it doesn’t now either, but I suppose it might have been messing with the parameter space.

As for looking at pairwise correlations between the covariates, there are few pairs of moderately highly correlated covariates, and as far as I can tell from looking at the marginal posteriors, what the model seems to do with each pair is to choose the stronger covariate & squeeze the other one towards 0. To me, that seems like a desirable behavior, at least for the purpose of building the lighter submodel with projection. Is that reasonable or is the collinearity something to worry about?

> [@avehtari](#):
>
> It’s difficult to see from the plot whether the full model performance is similar with RHS or normal prior?

Yes, based on LOO ELPD and RMSE (on the training data), the horseshoe performed about the same as the gaussian model does.

> [@avehtari](#):
>
> I think they are reasonably well behaved. In both cases there is similar bump in the curve, and it could bee good to check those predictor values.

The bump starts to appear after a somewhat unbalanced binary covariate (160 obs in one level, 433 in the other) is included (covariate #7). The other covariates before, and many after, are continuous. Could that be it?

---

<div class="post-metadata">

**Author:** ![abartonicek](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abartonicek/32/7020_2.png) [@abartonicek](https://discourse.mc-stan.org/u/abartonicek)\
**Post date:** [February 20, 2020, 10:14pm UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/11 "2020-02-20T22:14:54Z")

</div>

Thank you for replying so quickly!

---

<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, 2020, 9:29am UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/12 "2020-02-21T09:29:18Z")

</div>

> [@abartonicek](#):
>
> I’ve fixed it and now the selection trajectory looks like this:

Great.

> [@abartonicek](#):
>
> The unstandardized covariate did not appear until late in the feature selection trajectory, and it doesn’t now either, but I suppose it might have been messing with the parameter space.

Yes, it can.

> [@abartonicek](#):
>
> As for looking at pairwise correlations between the covariates, there are few pairs of moderately highly correlated covariates, and as far as I can tell from looking at the marginal posteriors, what the model seems to do with each pair is to choose the stronger covariate & squeeze the other one towards 0. To me, that seems like a desirable behavior, at least for the purpose of building the lighter submodel with projection. Is that reasonable or is the collinearity something to worry about?

projpred separates 1) making the best model with all covariates and 2) finding the smaller submodel with projection conditioanl on the best possible model from 1. When doing 1, we don’t need to try to squeeze many coefficient marginals towards 0 to make 2 to work well. When doing 1, we should use the best possible model and when we have correlating variables RHS is not exactly the correct model even if it can produce very good models. The reason why we use and recommend RHS even in the case of correlating covariates is that if there are only a small proportion of relevant covariates RHS is good for that and adding the extra part taking into account the correlations between those relevant covariates is complicated and based on the experiments has most of the time quite small effect on the model performance. If in your case Gaussian and RHS prior produce as good models in 1, then there shouldn’t be much difference in 2.

> [@abartonicek](#):
>
> Yes, based on LOO ELPD and RMSE (on the training data), the horseshoe performed about the same as the gaussian model does.

This was expected. You didn’t have very large number of covariates and you had n\>\>p, so you are not yet in the regime where you could more easily see differences between RHS and Gaussian prior.

> [@abartonicek](#):
>
> The bump starts to appear after a somewhat unbalanced binary covariate (160 obs in one level, 433 in the other) is included (covariate #7). The other covariates before, and many after, are continuous. Could that be it?

Probably. Adding a binary covariate can changesdistribution of individual predictive density values more and cross-validation can have higher variance for added binary covariates.

---

<div class="post-metadata">

**Author:** ![abartonicek](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/abartonicek/32/7020_2.png) [@abartonicek](https://discourse.mc-stan.org/u/abartonicek)\
**Post date:** [February 21, 2020, 11:36pm UTC](https://discourse.mc-stan.org/t/projpred-elpd-goes-down-and-rmse-goes-up-after-x-variables/13153/13 "2020-02-21T23:36:50Z")

</div>

Thank you so much for your time @avehtari! I’ve learned a lot, and I really appreciate the effort you go into in making awesome tools like `projpred` and helping other people to use them. Thank you @AlejandroCatalina as well for getting back to me about the error messages!