# LOO for Multivariate Probit

**URL:** <https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339>\
**Category:** Modeling\
**Tags:** loo\
**Created:** [February 14, 2022, 3:50pm UTC](https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339 "2022-02-14T15:50:17Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![fergusjchadwick](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fergusjchadwick/32/9765_2.png) [@fergusjchadwick](https://discourse.mc-stan.org/u/fergusjchadwick)\
**Post date:** [February 14, 2022, 3:50pm UTC](https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339/1 "2022-02-14T15:50:17Z")

</div>

Continuing the discussion from [Feedback request: Multivariate Probit Regression with GP](https://discourse.mc-stan.org/t/feedback-request-multivariate-probit-regression-with-gp/21517/7):

Brief problem summary: @bgoodri has developed this really handy implementation of the multivariate probit (code below). Unfortunately, it universally gets terrible k-hats in LOO. @martinmodrak made the excellent point in another thread and I wanted to continue discussion here to see if this is a viable way to get a robust LOO estimate (or whether this model is simply a poor candidate for LOO):

> Pareto-K is almost guaranteed to be bad if you just add up the contributions to `target` from this parametrization. In this case, you don’t compute the log likelihood of the observed values given the linear predictors. You actually compute the log likelihood given the linear predictors AND the nuisance parameters. Since the observed values have huge influence on the associated nuisance parameters, loo correctly treats them as having large influence on the model and thus having high k-hat.

Martin goes on to suggest:

> It could also make some kind of weird sense to compute the multivariate normal log likelihood of the nuisance parameters given the linear predictors and feed that to loo, but I can’t think completely clearly if that would correspond to a meaningful quantity or not.

I can very much see the logic behind this, but I don’t understand either LOO or the MVProbit implementation sufficiently to determine whether such a version of LOO would make sense (maybe @avehtari could comment?).

Thanks in advance!

Here is @bgoodri’s code (taken from [https://github.com/stan-dev/example-models/blob/master/misc/multivariate-probit/probit-multi-good.stan](https://github.com/stan-dev/example-models/blob/master/misc/multivariate-probit/probit-multi-good.stan)):

```no-highlight
data {
  int<lower=1> K;
  int<lower=1> D;
  int<lower=0> N;
  array[N, D] int<lower=0, upper=1> y;
  array[N] vector[K] x;
}
parameters {
  matrix[D, K] beta;
  cholesky_factor_corr[D] L_Omega;
  array[N, D] real<lower=0, upper=1> u; // nuisance that absorbs inequality constraints
}
model {
  L_Omega ~ lkj_corr_cholesky(4);
  to_vector(beta) ~ normal(0, 5);
  // implicit: u is iid standard uniform a priori
  {
    // likelihood
    for (n in 1 : N) {
      vector[D] mu;
      vector[D] z;
      real prev;
      mu = beta * x[n];
      prev = 0;
      for (d in 1 : D) {
        // Phi and inv_Phi may overflow and / or be numerically inaccurate
        real bound; // threshold at which utility = 0
        bound = Phi(-(mu[d] + prev) / L_Omega[d, d]);
        if (y[n, d] == 1) {
          real t;
          t = bound + (1 - bound) * u[n, d];
          z[d] = inv_Phi(t); // implies utility is positive
          target += log1m(bound); // Jacobian adjustment
        } else {
          real t;
          t = bound * u[n, d];
          z[d] = inv_Phi(t); // implies utility is negative
          target += log(bound); // Jacobian adjustment
        }
        if (d < D) {
          prev = L_Omega[d + 1, 1 : d] * head(z, d);
        }
        // Jacobian adjustments imply z is truncated standard normal
        // thus utility --- mu + L_Omega * z --- is truncated multivariate normal
      }
    }
  }
}
generated quantities {
  corr_matrix[D] Omega;
  Omega = multiply_lower_tri_self_transpose(L_Omega);
}

```

---

<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 16, 2022, 5:38pm UTC](https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339/2 "2022-02-16T17:38:23Z")

</div>

> [@fergusjchadwick](#):
>
> Unfortunately, it universally gets terrible k-hats in LOO. @martinmodrak made the excellent point in another thread

Yep, @martinmodrak’s point is accurate. You can see the same behavior when implementig overdispersed Poisson with n nuisance parameters [Roaches cross-validation demo](https://avehtari.github.io/modelselection/roaches.html#4_Poisson_model_with_%E2%80%9Crandom_effects%E2%80%9D)

I recommend cross-validating in the old fashioned way, that is, do K-fold-CV by re-fitting the model K times. There is vignette that might help [Holdout validation and K-fold cross-validation of Stan programs with the loo package • loo](https://mc-stan.org/loo/articles/loo2-elpd.html)

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [February 17, 2022, 12:11pm UTC](https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339/3 "2022-02-17T12:11:02Z")

</div>

> [@avehtari](#):
>
> I recommend cross-validating in the old fashioned way, that is, do K-fold-CV by re-fitting the model K times. There is vignette that might help [Holdout validation and K-fold cross-validation of Stan programs with the loo package • loo](https://mc-stan.org/loo/articles/loo2-elpd.html)

I tried reading how that’s done in the roaches example and I admit I am at loss on how the `kfold` method is able to compute the likelihood for the held-out observations - at some point the nuisance parameters (i.e. the per-row random effect) for the held-out observations have to be used, but it seems quite unclear to me how that happens… Are those sampled as new levels from the hyperprior? Is then a single sample used to compute log-likelihood? And should this be usable generally for all types of nuisance parameters?

---

<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, 2022, 5:42pm UTC](https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339/4 "2022-02-18T17:42:10Z")

</div>

> [@martinmodrak](#):
>
> I tried reading how that’s done in the roaches example and I admit I am at loss on how the `kfold` method is able to compute the likelihood for the held-out observations - at some point the nuisance parameters (i.e. the per-row random effect) for the held-out observations have to be used, but it seems quite unclear to me how that happens… Are those sampled as new levels from the hyperprior?

Yes. It’s part of rstanarm magic. rstanarm generates posterior draws for the random effect for a new group (this can be done in generated quantities by generating random draws from the population prior). rstanarm know that when it is predicting for new data with group factor that is not part of the original data, it uses the random draws for the new group. Able to use rstanarm predictions functions for new groups was there before kfold.

> [@martinmodrak](#):
>
> And should this be usable generally for all types of nuisance parameters?

If by “nuisance” you mean useful group level parameters (don’t hurt their feelings by calling them nuisance), then yes.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [February 20, 2022, 9:03pm UTC](https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339/5 "2022-02-20T21:03:58Z")

</div>

Thanks for the answer! So, to check my understanding, I’ll try to frame this in more general terms:

I have a model for data y\_1,...,y\_k, and two sets of parameters \theta and \nu\_1, ..., \nu\_k (all of which are potentially vectors). The likelihood for the full model decomposes as:

p(y\_1,...,y\_k | \theta) = \prod\_i p(y\_i | \theta) \\ p(y\_i | \theta) = \int p(y\_i | \theta, \nu\_i) p(\nu\_i | \theta) \mathrm{d}\nu\_i

Here, the \nu\_i are the definitely non-nuisance and very respecatble parameters.

When I am trying to perform k-fold CV, than after fitting without i-th part of the data I have posterior samples \theta\_{(-i),s}. I am interested in the quantity:

\log p(y\_i | y\_{-i}) = \log \int p(y\_i | \theta\_{-i}) p(y{-i} | \theta\_{-i}) \mathrm{d}\theta\_{-i} \approx \\ \approx \log \frac{1}{S} \sum\_s p(y\_i | \theta\_{(-i),s}) 

Where y\_{-i} is the data without i-th element.

If I understand Aki correctly, then I have two two options how to compute p(y\_i | \theta\_{(-i),s}):

**Variant 1 - Explicit integration:**

p(y\_i | \theta\_{(-i),s}) = \int p(y\_i | \theta\_{(-i),s}, \nu\_i) p(\nu\_i | \theta\_{(-i),s}) \mathrm{d}\nu\_i 

**Variant 2 - Sampling:**

1. For each s, draw a single new sample \nu\_{i,s} according to p(\nu\_i | \theta\_{(-i),s})
2. p(y\_i | \theta\_{(-i),s}) \approx p(y\_i | \theta\_{(-i),s}, \nu\_{i,s})

Obviously, when feasible, Variant 1 should always be 100% fine. If I understand Aki right, then either

a) Variant 1 and Variant 2 are different estimators of the same elpd, possibly differing only in variance of the estimates or  
b) (weaker variant) Variant 1 and Variant 2 will compute different values, but on average will lead to the same model comparison results.

Does that sound right? Am I missing something?

(I don’t need to see proofs or anything just trying to see if my understanding is correct)

Thanks very much!

> **A quick computational check supporting the answer a)**
>
> We’ll use the fact the neg. binomial can be rewritten as Poisson with gamma-distributed mean. Here is the Stan model `nb.stan`:
> 
> ```stan
> data {
> int<lower=0> N;
> int y[N];
> }
> 
> parameters {
> real<lower=0> mu;
> real<lower=0> phi;
> }
> 
> model {
> y ~ neg_binomial_2(mu, phi);
> }
> 
> generated quantities {
> vector[N] log_lik_explicit;
> vector[N] log_lik_sample;
> for(n in 1:N) {
> real gamma_var = gamma_rng(phi, phi);
> log_lik_sample[n] = poisson_lpmf(y[n] | mu * gamma_var);
> log_lik_explicit[n] = neg_binomial_2_lpmf(y[n] | mu, phi);
> }
> }
> 
> ```
> 
> Then we’ll compute 10-fold cross validation as well as directly using `loo` for a single fit:
> 
> ```r
> library(rstan)
> library(loo)
> mod_nb <- stan_model("nb.stan")
> 
> N <- 50
> y <- rnbinom(N, mu = 5, size = 1)
> 
> fold <- kfold_split_random(K = 10, N = N)
> 
> log_pd_kfold_explicit <- matrix(nrow = 4000, ncol = N)
> log_pd_kfold_sample <- matrix(nrow = 4000, ncol = N)
> 
> seed <- 1565233
> for(k in 1:10){
> data_train <- list(y = y[fold != k],
> N = sum(fold != k)
> )
> data_test <- list(y = y[fold == k], N = sum(fold == k))
> fit <- sampling(mod_nb, data = data_train, seed = seed, refresh = 0)
> gen_test <- gqs(mod_nb, draws = as.matrix(fit), data= data_test)
> log_pd_kfold_explicit[, fold == k] <- extract_log_lik(gen_test,parameter_name = "log_lik_explicit")
> log_pd_kfold_sample[, fold == k] <- extract_log_lik(gen_test,parameter_name = "log_lik_sample")
> }
> 
> (elpd_kfold_explicit <- elpd(log_pd_kfold_explicit))
> # Computed from 4000 by 50 log-likelihood matrix using the generic elpd function
> #
> # Estimate SE
> # elpd -124.1 7.6
> # ic 248.1 15.2
> 
> (elpd_kfold_sample <- elpd(log_pd_kfold_sample))
> # Computed from 4000 by 50 log-likelihood matrix using the generic elpd function
> #
> # Estimate SE
> # elpd -124.5 7.8
> # ic 248.9 15.6
> 
> # Comparison shows almost no difference
> (lc <- loo_compare(elpd_kfold_explicit, elpd_kfold_sample))
> # elpd_diff se_diff
> # model1 0.0 0.0   
> # model2 -0.4 0.3  
> 
> # Directly using loo
> fit_all <- sampling(mod_nb, data = list(y = y, N = N))
> # No problems with explicit form
> loo_explicit <- loo(fit_all, pars = "log_lik_explicit")
> # All k-hats are crazy high when using samples
> loo_sample <- loo(fit_all, pars = "log_lik_sample")
> 
> ```

---

<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, 2022, 5:03pm UTC](https://discourse.mc-stan.org/t/loo-for-multivariate-probit/26339/6 "2022-02-21T17:03:49Z")

</div>

> [@martinmodrak](#):
>
> Variant 1 - Explicit integration:

Yes, see, e.g.

- [Aki Vehtari, Tommi Mononen, Ville Tolvanen, Tuomas Sivula and Ole Winther (2016). Bayesian leave-one-out cross-validation approximations for Gaussian latent variable models. _Journal of Machine Learning Research_ , 17(103):1−38.](http://jmlr.org/papers/v17/14-540.html)

> [@martinmodrak](#):
>
> Obviously, when feasible, Variant 1 should always be 100% fine.

Define fine?

> [@martinmodrak](#):
>
> a) Variant 1 and Variant 2 are different estimators of the same elpd, possibly differing only in variance of the estimates or

Yes, when you do leave-one-out CV.

It’s good to make difference to two variants discussed in

- [Merkel, Furr, and Rabe-Hesketh (2019). Bayesian Comparison of Latent Variable Models: Conditional Versus Marginal Likelihoods. _Psychometrika_ 84:802-829.](https://link.springer.com/article/10.1007/s11336-019-09679-0)

where the other variant corresponds to leave-one-group-out which is different if each latent parameter has more than one observation.